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HIGHLIGHTS 


►  A  theory  of  coupled  mechano-diffusional  driving  forces  for  fracture  (/-integral)  is  developed. 

►  Analysis  shows  that  f  and  K  are  not  uniquely  related  due  to  full  deformation-diffusion  coupling. 

►  Global  yielding  and  lithiation-induced  softening  reduce  energy  release  rate  in  thin-film  Li/Si  electrodes. 

►  Operation  at  higher  Li  concentrations  mitigates  fracture  of  thin-film  Li/Si  electrodes. 

►  A  design  map  is  developed  for  guiding  the  safe  operation  of  thin-film  Li/Si  electrodes. 


ARTICLE 


N  F  O 


A  B  S  T  R 


C  T 


Article  history: 

Received  31  October  2012 
Received  in  revised  form 
28  November  2012 
Accepted  6  December  2012 
Available  online  13  December  2012 


Keywords: 

Plasticity 

Energy  release  rate 
Stress  intensity  factor 


Li/Si  undergoes  significant  softening  in  the  form  of  rapid  decreases  in  elastic  modulus  and  yield  stress  as 
lithium  concentration  increases.  To  investigate  how  this  lithiation-induced  softening  affects  the  fracture 
behavior  of  electrodes,  we  formulate  a  /-integral  for  coupled  mechanical  deformation  and  mass  diffusion 
processes.  This  measure  is  used  to  analyze  mechano-diffusional  driving  forces  for  fracture  through 
simulations  using  a  mixed  finite  element  framework.  Calculations  show  that  under  tensile  loading,  Li 
accumulates  in  front  of  crack  tips,  leading  to  an  anti-shielding  effect  on  the  energy  release  rate.  For 
a  pre-cracked  Li/Si  thin-film  electrode,  it  is  found  that  the  driving  force  for  fracture  is  significantly  lower 
when  the  electrode  is  operated  at  higher  Li  concentrations  —  a  result  of  more  effective  stress  relaxation 
via  global  yielding.  The  results  indicate  that  operation  at  higher  concentrations  is  an  effective  means  to 
minimize  failure  of  thin-film  Li/Si  alloy  electrodes.  A  design  map  for  avoiding  failure  is  developed. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Electrode  degradation  due  to  internal  stress  build-up  during 
cycling  has  been  one  of  the  key  challenges  for  secondary  battery 
designers.  When  guest  atoms  are  inserted  or  extracted  from  the 
host  in  an  electrode,  the  material  in  the  electrode  expands  or 
contracts,  inducing  stresses  that  may  cause  material  cracking  [1], 
This  issue  is  especially  important  in  alloy-based  electrode  mate¬ 
rials,  which  are  attractive  because  their  capacities  are  much  higher 
than  that  of  graphite,  but  at  the  same  time  experience  much  higher 
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volume  expansions.  Size  reduction  has  been  proven  to  be  an 
effective  means  to  mitigate  mechanical  failure  [2],  In  particular, 
alloy-based  electrodes  made  of  nano-sized  structures  have  been 
demonstrated  to  significantly  improve  cyclability  [3—7], 

Many  models  have  been  proposed  to  characterize  the  buildup 
and  mitigation  of  stresses  in  Li-ion  battery  electrodes  [2,8-14], 
Bower  et  al.  [15]  developed  a  comprehensive  framework  and  used 
it  to  analyze  time-dependent  plasticity  in  thin-film  Li/Si.  The  effect 
of  stress-enhanced  diffusion  (SED)  was  analyzed,  revealing  signif¬ 
icant  reduction  in  stress  due  to  a  mechanical  driving  force  for 
diffusion  when  the  deformation  is  in  the  elastic  regime  [16],  This 
enhancement  effect  can  nevertheless  be  diminished  or  even  fully 
reversed  when  the  material  deforms  plastically  [17],  Deshpande 
et  al.  [18]  considered  the  effect  of  surface  stresses  and  concluded 
that  surface  effect  can  reduce  the  tensile  stresses  in  nano-sized 
electrodes,  thereby  improving  electrode  cyclability.  Zhao  et  al. 
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[19]  considered  plastic  deformation  and  showed  that  inelastic  flow 
can  significantly  alleviate  stresses  in  Li/Si.  The  mechanisms  for 
stress  reduction  considered  in  these  analyses  include  stress- 
enhanced  diffusion  of  Li,  surface-effect-induced  compressive 
stresses,  and  plasticity. 

Among  the  many  candidates  for  next-generation  negative 
electrode  materials,  silicon  is  considered  one  of  the  most  promising 
because  of  its  high  capacity  and  low  cost.  Although  pure  silicon  is 
brittle,  both  experimental  and  theoretical  studies  show  that  lithi- 
ated  silicon  can  undergo  significant  plastic  deformation  [20—22], 
Studies  also  show  that  such  a  transition  from  brittle  to  ductile  is 
gradual,  with  both  the  elastic  modulus  and  yield  stress  decreasing 
continuously  as  Li  concentration  increases  [20,23,24],  However, 
how  such  lithiation-induced  softening  affects  the  fracture  tendency 
of  Li/Si  electrodes  is  neither  well-understood  nor  quantified. 

The  fact  that  lithiated  Li/Si  deforms  plastically  indicates  that 
elastic-plastic  fracture  mechanics,  instead  of  linear  elastic  fracture 
mechanics  (LEFM),  should  be  used  to  characterize  failure.  What 
complicates  the  problem  is  the  fact  that  the  distribution  of  Li,  which 
gives  rise  to  stresses,  is  in  turn  affected  by  the  stress  state 
[15,25,26],  When  a  crack  is  mechanically  loaded,  Li  ions  would 
accumulate  in  front  of  the  crack  tip,  in  a  fashion  similar  to  that  seen 
hydrogen  embrittlement  problems  [27,28],  This  accumulation  of  Li 
ions  at  the  crack  tip  can  have  profound  implications  on  the  driving 
forces  for  fracture. 

Recently,  Ryu  et  al.  [11]  proposed  a  framework  for  calculating 
the  energy  release  rate  /  for  cracks  in  Si  nanowire  electrodes  and 
used  the  framework  to  study  the  size  dependence  of  fracture.  This 
theory  relies  on  an  “effective  diffusivity”  Deg  into  which  the  effect 
of  stress  on  diffusion  is  lumped.  Such  an  effective  diffusivity  is  only 
applicable  to  linear  elastic  problems.  Specifically,  when  a  material 
deforms  plastically,  diffusion  is  no  longer  governed  by  DeJy  and  the 
overall  diffusion  kinetics  depends  on  the  nature  of  mechanical 
constraints  as  well  as  other  factors  [17],  Even  in  the  elastic  regime, 
calculations  based  on  Deg  can  only  capture  stress-induced 
enhancement  to  diffusion  [16]  in  regions  far  from  the  crack  tip 
but  cannot  capture  the  trapping  of  Li  ions  at  the  crack  tip  [28].  To 
correctly  account  for  the  underlying  processes  responsible  for 
fracture  which  include  Li  trapping  and  plasticity,  a  fully-coupled 
theory  is  needed  to  assess  the  energy  release  rate. 

Here,  such  a  fully-coupled  theory  is  developed  and  used  to 
analyze  the  coupled  mechano-diffusional  driving  forces  for  fracture 
in  electrode  materials.  The  analyses  take  advantage  of  a  mixed 
finite  element  framework  [29,30],  We  first  formulate  a  /-integral  for 
coupled  mechanical  deformation  and  mass  diffusion  processes  as 
a  measure  for  the  driving  forces  for  crack  growth.  By  treating  the 
crack  tip  as  a  separate  thermodynamic  system,  the  formulation 
entails  detailed  account  of  balances  of  mass  and  energy  and 
evolution  of  entropy.  All  relevant  energy  forms  that  may  contribute 
to  fracture  are  explicitly  tracked.  It  is  found  that  the  standard  form 
of /-integral  for  energy  release  is  no  longer  path-independent  when 
coupled  mechano-diffusion  driving  forces  are  present.  Instead,  an 
area  integral  similar  to  those  in  hygrothermal  [31]  and  dynamic 
[32—34]  problems  must  be  included.  A  numerical  scheme  for 
implementing  this  path-independent  formulation  through  finite 
element  simulations  is  also  given. 

The  numerical  study  based  on  the  framework  is  carried  out  in 
a  progressive  manner,  with  increasing  complexity  in  each  step.  In 
Section  3.1,  we  first  consider  the  linear  elastic  case  with 
concentration-independent  properties,  so  that  the  near-tip  stress 
fields  and  the  energy  release  rate  can  be  compared  with  those  given 
by  purely  mechanical  LEFM.  It  is  found  that  stress-induced  lithium 
redistribution  significantly  affects  energy  release,  but  has  no  effect 
on  in-plane  stress  fields  and  hence  the  stress  intensity  factor  Kj.  The 
implication  of  such  a  dichotomy  is  discussed  by  drawing  an  analogy 


to  the  difference  between  plane-strain  and  plane-stress  situations 
in  LEFM. 

The  next  step  in  the  analyses  (Section  3.2.1 )  assumes  a  finite  yield 
stress  so  that  the  material  can  deform  plastically.  This  scenario 
involves  elasto-plastic  deformation  with  concentration-independent 
elastic  properties  and  yield  stress.  The  calculation  focuses  on 
a  surface  crack  in  a  thin-film  electrode.  Depending  on  whether  the 
electrode  yields  globally,  two  regimes  of  response  can  be  identified. 
Before  global  yielding,  the  energy  release  rate  is  governed  by  an 
effective  crack  length  which  is  analogous  to  that  proposed  by  Irwin. 
After  the  global  yielding  threshold,  however,  it  is  found  that  plastic 
deformation  significantly  reduces  the  energy  release  rate  and  softer 
materials  generally  have  much  lower  energy  release  rate  levels. 

Finally  in  Section  3.2.2,  the  effect  of  lithiation-induced-softening 
on  fracture  tendency  is  added.  The  calculations  consider 
deformation/diffusion  coupling,  plastic  flow,  and  lithiation-induced 
material  softening  through  decreases  in  the  elastic  modulus  and 
yield  stress  as  Li  concentration  increases.  The  analyses  lead  to  a  design 
map  for  the  configuration  and  operation  of  thin-film  Li-Si  electrodes. 

2.  Energy  release  rate  under  full  diffusion-deformation 
coupling 

We  consider  the  energy  release  rate  for  fracture  in  a  solid 
electrode  with  full  diffusion— deformation  coupling.  The  material, 
which  is  composed  of  two  chemical  species:  host  (denoted  as  H) 
and  guest  (lithium),  is  assumed  to  be  highly  conductive  so  that  the 
whole  electrode  is  at  the  same  electric  potential.  The  electrode  is 
assumed  to  be  fully  amorphized  from  the  beginning  so  that  phase 
separation  due  to  crystalline/amorphous  transition  does  not  need 
to  be  considered.  During  charge/discharge,  lithium  atoms  diffuse, 
while  the  host  atoms  are  assumed  to  undergo  convection  but  not 
diffusion  [29],  The  movement  of  the  host  atoms  therefore  defines 
the  continuum  deformation  x  =  x(X,  t).  The  Eulerian  concentra¬ 
tions  of  the  host  and  lithium,  namely  the  atomic  numbers  per  unit 
current  volume,  are  denoted  by  cH  and  cLl,  respectively.  Their 
Lagrangian  counter  parts  CH  and  CLi  are  related  to  cH  and  cLi 
through  CH  =  det(F )cH  and  CLi  =  det(F )cLi,  where  F  =  3x/3X  is 
the  deformation  gradient.  Under  the  assumption  of  zero  host 
diffusivity,  CH  does  not  change  with  time,  and  a  single  dimen¬ 
sionless  composition  ? =CU/CH  can  be  used  to  conveniently  char¬ 
acterize  the  local  state  of  charge.  This  dimensionless  composition 
?  corresponds  to  the  lithium  number  per  host  in  the  chemical 
formula  LifH.  The  maximum  ?  at  the  fully  charged  state  is  denoted 
as  ?max.  and  the  state  of  charge  (SOC)  relative  to  this  fully  charged 
state  is  denoted  as  2=f/?max. 

Following  the  standard  theory  of  large-deformation  plasticity, 
a  Lee-type  decomposition  [15,35]  can  be  performed  for  the  defor¬ 
mation  gradient,  i.e., 

F  =  Fe-Fsp-F p,  (1) 

where  F6,  Fsr  and  Fp  are  the  elastic,  stress-free  volumetric  and 
plastic  deformation  gradients,  respectively  [12,15],  The  associated 
rates  of  deformation  are 


DP=  [Fe-FSF.Fp-(FP)-1-(Fsf)  1  -(F6)-1]  m,and 
D  =  [P(F)_1]  sjmj  =  De  +  Dsf  +  DP. 
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We  assume  the  stress-free  expansion  due  to  concentration 
change  is  isotropic  so  that 

Fsf  =  VSF(?)I  =  [jsf(?)]1/3I,  (3) 

where  JSF(?)  =  det(FSF)  measures  the  volume  change  due  to 
changes  in  Li  concentration.  V^Esf®^)]1/3  is  the  associated 
stretch. 

The  diffusion— deformation  coupling  comprises  of  two  interde¬ 
pendent  aspects.  The  first  aspect,  the  effect  of  diffusion-induced 
stress  (DIS),  arises  due  to  a  C  -dependent  term  in  the  constitu¬ 
tive  relationship  between  D  and  a  [cf.  Eqs.  (A.4)  and  (A.5)],  where 
C  is  the  time  rate  of  CLl  and  <r  is  an  objective  rate  (e.g.,  the  Jau- 
mann  rate)  of  the  Cauchy  stress  cr.  The  second  aspect,  the  effect  of 
mechanical  driving  forces  on  diffusion,  is  embodied  by  a  stress- 
driven  term  ]L,-mech  jn  the  total  Lagrangian  flux  ]h  of  lithium  [cf. 
Eq.  (A.17)].  These  two  aspects,  together  with  the  conservation  laws 
of  mass  [Eq.  (A.9)]  and  momentum  [Eq.  (A.7)],  govern  the  material 
response  in  a  solid  electrode  during  charge  and  discharge.  A  brief 
overview  of  this  continuum  framework  is  given  in  Appendix  A. 

The  consideration  of  the  energy  release  rate  (/)  for  fracture  is 
based  on  an  account  of  the  balance  of  mechanical  work  in  a  cracked 
body  [33].  For  electrodes  whose  response  is  affected  by  diffusion- 
deformation  coupling,  care  must  be  taken  because  mechanical 
energy  is  no  longer  the  only  energy  source  at  work.  Specifically, 
chemical  energy  which  drives  the  mixing  between  Li  and  host  may 
also  constitute  a  significant  contribution  to  the  overall  energy 
balance.  Such  coupled  driving  forces  for  fracture  in  electrode 
materials  are  similar  to  those  in  hygro-thermal  systems.  Chen  [31  ] 
proposed  a  generalized  contour  integral  method  to  evaluate  the 
latter.  The  method,  however,  is  path-dependent  and  requires 
infinitesimal  contours  near  the  crack  tip  for  accurate  evaluation. 

In  this  section,  we  use  non-equilibrium  thermodynamics  to 
develop  a  path-independent  /-integral  for  the  coupled  mechano- 
diffusional  driving  forces  for  fracture  in  electrode  materials  that 
undergo  large  elasto-plastic  deformations.  The  analysis  treats  the 
crack  tip  as  a  separate  thermodynamic  system,  therefore  allowing 
a  systematic  account  of  the  balance  of  mass  and  energy  and  the 
production  of  entropy.  A  numerical  scheme  for  implementing  this 
path-independent  formulation  in  finite  element  simulations  is  also 
developed. 

2.1.  Energy  balance  and  entropy  production 

The  analysis  of  energy  release  rate  under  full  diffusion- 
deformation  coupling  is  based  on  accounts  of  energy  and  entropy. 
For  any  material  point  in  the  continuum  body,  the  first  law  of 
thermodynamics  [36]  requires  that 

e  =  Vx-(<rFF1-v-f-uFjFi),  (4) 

where  e  is  the  internal  energy  per  unit  reference  volume, 
v  =  dx/dt|x  is  particle  velocity,  aPK1  =  det(F)F_1  -  a  is  the  first 
Piola-Kirchhoff  stress,  and  uLi=de/dCLi\e  is  the  partial  atomic 
energy  of  Li.  ]q  and  Ju  are  the  heat  flux  and  Li  atom  flux,  respec¬ 
tively,  both  measured  in  the  reference  configuration.  In  writing  Eq. 
(4),  we  have  assumed  that  there  is  no  external  volumetric  heat 
source.  Under  this  condition,  the  2nd  law  of  thermodynamics  can 
be  stated  as  [36] 


Here,  rj  and  p  are  the  Lagrangian  density  and  flux  of  entropy, 
respectively;  8  is  temperature;  and  2  is  the  rate  of  entropy 
production  per  unit  reference  volume.  In  this  paper,  isothermal 
conditions  are  assumed  to  prevail  [i.e.,  8  =  6(X.  t)  =  constant], 
therefore 

2  =  V  +  ivx.  [j*  +  (uLi  -  mu)J°]  >  0.  (6) 

The  entropy  production  rate  2  phenomenologically  character¬ 
izes  the  irreversible  nature  of  the  processes  in  the  continuum.  For 
materials  that  deform  plastically,  2  can  be  identified  as  the  sum  of 
2P  (which  is  due  to  plastic  flow)  and  2NP  (which  is  due  to  dissipa¬ 
tion  mechanisms  other  than  plastic  deformation,  e.g.  diffusion),  i.e., 

2  =  2P  +  2np, 

2P  =  .  dp  >  0,  and  1  (7) 

2np  >  0.  | 

2.2.  The  crack  tip  subsystem 

Following  Moran  and  Shih  [37,38]  and  Freund  [39],  we  consider 
a  2-D  body  in  the  (ei,  e2)  plane  with  a  crack  extending  along  the 
e,  direction,  as  illustrated  in  Fig.  la.  At  time  t,  the  crack  has  length 
l  =  l(t )  in  the  Lagrangian  configuration.  The  crack  tip,  which 
propagates  at  speed  1  in  the  reference  configuration,  can  be  isolated 
using  a  small  Lagrangian  contour  T  which  translates  along  with  the 
crack  tip.  The  nominal  out-of-plane  thickness  (i.e.,  nominal  thick¬ 
ness  along  e3)  of  the  2-D  body,  also  measured  in  the  reference 
configuration,  is  denoted  as  dcrack. 

Consider  a  larger  contour  Co  which  is  fixed  with  respect  to  the 
material.  There  is  no  singularity  in  domain  A  bounded  by  T,  C0  and 
the  two  crack  surfaces  C_  and  C+.  In  the  crack-tip  domain  Vr 
bounded  by  T,  C_  and  C+,  however,  singularity  exists.  We  treat  this 
singular  point  in  V,  as  a  stand-alone  thermodynamic  subsystem  y 
which  can  exchange  mass,  energy  and  entropy  with  the 
surrounding  continuum  body.  The  Helmholtz  free  energy  of  this 


2  =  i)  +  VX\T  >  0,  where 


Fig.  1.  Energy  balance  at  the  crack  tip,  reckoned  in  the  Lagrangian  configuration,  (a) 
Domain  A  is  enclosed  by  T  and  C  =  C_  +  Qj  +  C+.  Nj  and  Mj  are  unit  vectors  normal  to 
(5)  T  and  Co,  respectively.  Vr  is  the  domain  enclosed  by  T,  C_  and  C+.  The  infinitesimal 
region  Vr  exchanges  energy  and  entropy  with  the  crack  subsystem  which  is  charac¬ 
terized  by  its  free  energy  (b)  A  special  type  of  contour  r  with  rectangular  shape. 
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subsystem  per  unit  out-of-plane  thickness  dcrack  is  assumed  to  take 
the  form  of 

0crack  =  Qcrack^jcrack  ^  (8) 


Vr.  In  the  derivation  of  Eq.  (12),  use  has  been  made  of  the  Reynold 
transport  theorem  on  the  local  forms  of  Eqs.  (A.9),  (4)  and  (5). 

Combining  the  last  two  equations  of  Eq.  (12)  with  Eq.  (9)  under 
isothermal  conditions  gives 


where  tfmck  is  the  temperature  and  nu  is  the  amount  of  Li  per  unit 
dcmck.  It  turns  out  that  Eq.  (8)  completely  characterizes  the  ther¬ 
modynamic  properties  of  the  crack  subsystem.  Specifically,  the 
surface  energy  per  unit  crack  area  (y),  the  entropy  per  unit  out-of¬ 
plane  thickness  {Scrack),  the  chemical  potential  for  Li  (jxL,  crack)  and 
the  internal  energy  per  unit  out-of-plane  thickness  (Ucrack)  can  all 
be  uniquely  determined  via  <Pcrack  as 


1  ;)<j)crack 

T  “  2  dl 


d([>crack 

~d6crack’ 


Li, crack 


d<PCrack 
dnLi, crack’ 


and 


Ucrack  =  <Pcmck  +  scrack8crack, 


(9) 


respectively.  Since  isothermal  conditions  are  assumed,  temperature 
is  uniform  in  the  continuum  domain  and  the  crack  subsystem  y  is  in 
thermal  equilibrium  with  its  surrounding  material  in  Vr,  i.e. 


f)crack  =  d{X,  t)  =  const.  (10) 

Under  this  assumption,  the  incremental  form  of  Eq.  (8)  can  be 
simplified  with  the  help  of  Eq.  (9)  into 

0CmCk  =  2y(  +  n^i, crack  j^Li, crack  (U) 

2.3.  Energy  release  rate 

Without  loss  of  generality,  we  assume  that  y  can  exchange  mass, 
energy  and  entropy  with  the  continuum  body  (domain  A)  only 
through  Vr.  In  other  words,  we  neglect  any  direct  mass,  energy  and 
entropy  exchange  between  y  and  A  (one  possible  mechanism  for 
such  direct  exchange  is  surface  adsorption  on  C_  and  C+).  Under 
this  assumption,  the  balances  of  mass,  energy  and  entropy  for  the 
combined  system  of  y  and  Vr  (denoted  as  y  +  Vr )  require  that 


d>+!  J  4>dA0  =  J of\Njdr0  +  J WyNjdTo 

Vrr  r  r  Tr  1  (13) 

-  /  nLiJLi -Mr o  -  8  2crack  +  /  2dA0  , 

r  Vr 

where  0  =  c-r]8  is  the  Helmholtz  free  energy  of  the  continuum 
per  unit  reference  volume.  Again,  scrack  >  0  and  2  >  0  here, 
according  to  the  2nd  law  of  thermodynamics. 

Equation  (13)  embodies  the  balance  that  accounts  for  the 
transfer  and  dissipation  of  Helmholtz  free  energy  around  the  crack 
tip.  The  significance  of  its  terms  can  be  explained  as  follows.  The 
left  hand  side,  <t>  +  ( d/dt )  fv  </>dA0,  is  the  time  rate  of  the  total 
Helmholtz  free  energy  in  the  combined  system  y  +  Vr.  This  system 
comprises  of  everything  enclosed  by  the  moving  boundary 
T,  including  the  singularity  point.  The  four  terms  on  the  right  hand 
side  are,  respectively,  the  rates  of  the  mechanical  work  done  by 
domain  A  to  y  +  Vr  (1st  term),  the  free  energy  swept  into  y  +  Vr  by 
the  moving  boundary  T  (2nd  term),  the  free  energy  conveyed  into 
y  +  Vr  by  mass  flux  (3rd  term),  and  the  loss  of  available  free  energy 
in  y  +  Vr  due  to  irreversible  dissipation  (4th  term).  Here,  the  shape 
of  T  is  arbitrary  and,  except  Eq.  (8),  no  specific  assumptions 
on  material  constitutive  behavior  are  made  in  the  derivation  of 
Eq.  (13). 

Note  that  Eq.  (13)  accounts  for  different  dissipation  mechanisms 
in  the  continuum,  including  inelastic  flow  [e.g„  2P  in  Eq.  (7)]  and 
mass  transport  [which  can  be  lumped  into  2NP  in  Eq.  (7)].  The  only 
requirement  is  that  the  2nd  law  of  thermodynamics  is  satisfied 
such  that  2crack  >  0  and  2  >  0.  Now,  consider  the  specific  consti¬ 
tutive  response  described  in  Appendix  A,  in  the  limit  that  T  shrinks 
to  the  crack  tip.  For  the  form  of  0  given  by  Eqs.  (A.12)-(A.15),  the 
singularity  satisfies  the  condition 

ftl(af/*Mo)  ‘  0  (,4) 


mass :  hLUmck  +  ^  j CLidA0  =  -JffNjdTo  +  j CLi'l8vNjdr0, 
vT  r  r 

energy  :  Ucmck  +~  J  edA0  =  J  [ojf1  vt -jf  -  uLiJjU  +  eidy  .NjdF0: 

Vr  T 

entropy:  Sarack  +  ^  J  r]dA0  =  2crack  +  J  2dA0  ~  J  +  ^  g  ^  ^1+vl5\j^jNJdr0. 


Here,  /v  (•')dAo  stands  for  area  integral  over  the  2D  Lagrangian 
domain  Vr,  while  [r(»)dr0  stands  for  path  integral  along  the  2D 
contour  T.  The  area  element  dA0  and  line  segment  dF0  are  both 
measured  in  the  reference  configuration.  zcrack  >  0  is  the  entropy 
production  rate  (per  unit  d0'11^),  which  can  be  due  to  dissipation 
inside  y  or  dissipation  associated  with  interactions  between  y  and 


This  is  true  as  long  as  a  <  1  for  the  asymptotic  behavior  of 
(jy  ~  1  /ra  when  r  -►  0,  where  r  is  the  radial  distance  measured  from 
the  crack  tip.  On  the  other  hand,  the  form  of  chemical  potential  gLi 
used  in  Eq.  (A.13)— (A.15)  ensures  that  0  <  ?  <  ?max  even  when 
dm  ~  l/r“  goes  to  infinity  as  r— >0.  The  fact  that  ?  (and  hence  CLl)  is 
continuous  and  bounded  leads  to 
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f  ,'Tofs  JcLidAo)  =  °’and 

J  V  Vr  7 
\im(  J  CLil5vNjdr0\  =  0. 


Mass  conservation  near  the  crack  tip  [cf.  Eq.  (12)]  when  r ->0, 
therefore,  takes  the  form 


hLl’ crack  _  -  ]irr^  jjj'NjdI'0,  (16) 

r 

i.e.,  Li  atoms  that  transport  into  y  +  Vr  through  F  are  predomi¬ 
nantly  stored  in  the  crack  tip  subsystem  y  if  f->0. 

Generally  speaking,  lithium  chemical  potential  gLi  in  Vr  is 
position-dependent  and  not  necessarily  equal  to  jxLi  crack  in  the  y 
subsystem.  When  pu=£pLi’crack,  the  transport  of  Li  into  y  (at  the  rate 
of  fiu,crack )  is  dissipative  and  can  be  characterized  by  a  positive- 
definite  dissipation  rate  of  JimJnLLCrack(/j,L'  —  gIi’cr “*)],  For 
simplicity,  such  an  dissipation  mechanism  is  not  considered  in  this 
paper,  i.e.,  it  is  assumed  that  y  is  in  chemical  equilibrium  with  its 
immediate  surroundings  such  that 


depend  on  the  state  of  charge  (i.e.,  Li  concentration)  in  the  elec¬ 
trode.  Linder  this  assumption,  the  critical  condition  for  isothermal 
quasi-static  crack  growth  accounting  for  deformation-diffusion 
coupling  can  be  stated  as 

J  =  Jcr  >  2y.  (22) 

Equation  (22)  is  the  Griffith  condition  for  fracture  growth  when 
diffusion-deformation  two-way  coupling  is  present.  It  should  be 
noted  that  as  a  fracture  mechanics  criterion,  the  Griffith  condition 
generally  does  not  apply  to  fatigue  crack  growth  and  the  failure  of 
secondary  battery  electrodes  can  occur  under  static  or  cyclic  loads 
(charge-discharge).  However,  for  an  electrode  to  last,  it  has  to 
survive  the  first  few  cycles  of  loading  during  which  fracture  is 
governed  by  Eq.  (22).  As  will  be  shown  in  the  numerical  analysis  in 
Section  3,  even  this  requirement  alone  puts  significant  constraints 
on  the  design  of  battery  electrodes.  We,  therefore,  use  Eq.  (22)  as 
the  criterion  for  battery  failure  in  this  paper,  and  leave  the 
discussion  of  cyclic  fatigue  to  future  studies.  It  should  be  also  noted 
that  even  for  situations  in  which  fatigue  is  important,  analyses  such 
as  the  Paris-law  generally  require  the  calculation  of  fracture  driving 
force  (in  the  form  of  K )  or  J).  The  fully-coupled  theory  in  this  paper 
can,  therefore,  also  be  regarded  as  an  essential  part  of  future  models 
for  cyclic  failure  in  battery  electrodes. 


^Li, crack  =  lim  ^ Li 


(17) 


2.4.  Path-independent  integral  for  energy  release  rate 


This  assumption  [Eq.  (17)]  and  Eq.  (16)  immediately  lead  to 

fint  J pLif  •  Nd/ o  p;  ljiLi’crackfiLi’crack.  (18) 

r 

Combining  Eqs.  (11),  (13),  (14)  and  (18),  therefore,  yields 


The  specific  form  of  Eq.  (20)  is  general  but  inconvenient  for  the 
numerical  evaluation  of  energy  release  rate  because  the  integral 

/(-T) «  J  [M,  -  NjdF0  (23) 

r  L 


firn  J  Vi  +  </d(5qj  NjdF0  =  2y  1  +  6  \^scrack  +  firn  J SdA0J . 

r  Mr 

(19) 

The  left  hand  side  of  Eq.  (19)  is  the  time  rate  at  which  energy  is 
imparted  into  the  crack  tip  region  during  crack  growth.  Under  the 

conditions  of  local  steady  state,  Vj - i(dUi/dX-,).  The  driving  force 

for  crack  growth,  namely  energy  provided  per  unit  new  crack  area 
created,  is  therefore  given  by 

JB  rP0  J / +  <Pi5v]N]dI'o 

/r  a,  i  (2°) 

-S/b-“T5 x)]NA°- 

The  right  hand  side  of  Eq.  (19),  on  the  other  hand,  comprises  two 
parts.  The  first  part  2yi  corresponds  to  the  increase  of  free  energy  in 
the  subsystem  y,  while  the  second  part  d[Xcmck  +  Hrn  /V(  i’cM0] 
characterizes  the  free  energy  loss  due  to  dissipation.  As  for  any 
irreversible  processes,  this  dissipation  is  generally  history- 
dependent,  in  the  sense  that  the  combined  resistance 

e^Ecrack+  lim  JzdAo^ 

Jcr  =  27  +  Hm - ;  V' -  (21 ) 

depends  on  the  history  of  crack  growth  even  for  the  same  ther¬ 
modynamic  state  (1, 0cracfc,  nLl)  of  the  crack  tip  subsystem  y.  For 
simplicity,  the  standard  assumption  of  fracture  mechanics  that 
Jcr  —  Jcr(I)  is  used  here,  with  the  understanding  that  Jcr  may  also 


is  path-dependent  for  coupled  diffusion/deformation  problems. 
Indeed,  for  any  continuum  domain  Q  bound  by  dQ,  it  can  be  proven 
that  (see  Appendix  B) 

J*(dQ)  =  -  J wp5yNjdF0  +  J PLi^dA0,  (24) 

afi  a 


where  the  plastic  potential  vi/1  is  defined  as 


Here,  aPK2  =  det(F)F_1 -ct-F_t  is  the  2nd  Piola-Kirchhoff 
stress,  and  Ep  is  the  Lagrangian  strain  of  plastic  deformation.  Sub¬ 
tracting  Eq.  (24)  from  Eq.  (23)  leads  to  the  path-independent 
integral  form  of 

J(r)  =  J  [4>0y  ~  tfPidUt/dX^Njdro 


The  integral  J(T)  is  path-independent,  in  the  sense  that 
J(dfl)  =  0  for  the  boundary  dQ  of  any  continuum  domain  Q.  On  the 
other  hand,  contour  T  in  Eq.  (20)  can  be  of  arbitrary  shape  as  long  as 
its  size  is  infinitesimal.  To  establish  a  link  between  Eqs.  (20)  and  (26), 
consider  the  special  type  of  T  shown  in  Fig.  1(b).  As  pointed  out  by 
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Freund  [39],  if  the  rectangular  contour  is  shrunk  onto  the  crack  tip  by 
first  letting  $2  — >0  and  then  <5 1  ->0,  one  has 

'■RR/ b-^'S^Ao.  (27) 

Such  rectangular  contours  are  especially  convenient  because 
horizontal  segments  2  and  4  do  not  contribute  to  the  integral 
Ji_2_3-4-5WP<5i/Njdr0  [cf.  Fig.  2(b)].  Even  if  wp  is  singular  at  the 
crack  tip,  the  plasticity  term  in  Eq.  (26)  still  vanishes  as  long  as 
appropriate  limiting  process  is  considered,  i.e. 

lim  lim  [ wp5yNjdr0  =  0.  (28) 

a,— o  (52->o  7 
r 

On  the  other  hand,  for  domain  Vr  bound  by  T  [cf.  Fig.  2(b)], 
/v  fiLt(dCLl/dXj)dAocan  be  evaluated  and  goes  to  0  when  <52  — >  0  and 
<5/->0,  i.e., 

Urn  hrn  J ftLi(dCLi/dX^dA0  =  0.  (29) 

Combining  Eqs.  (26)— (29),  therefore,  leads  to  the  path- 
independent  ./-integral  for  coupled  mechano-diffusional  driving 
forces  in  the  form  of 

I=m  =  /[(0  +  wp)5u-of1M]Njdro  -  / ^0, 

T  WT 

(30) 

where  T  is  arbitrary.  With  Eqs.  (A.11)-(A.13)  and  some  algebra,  Eq. 
(30)  can  be  recast  into 

J  =m  =  J  [  (w  +  wP)  5y  -  (Tf1  (dUi/dXj )]  Njdr0 
-/($! 

where  w  is  the  elastic  strain  energy  given  by  Eq.  (A.12). 


Equations  (30)  and  (31)  will  be  used  to  evaluate  energy  release 
rates  from  numerical  results  given  by  finite  element  simulations. 
The  main  advantage  of  this  path-independent  form  lies  in  the  fact 
that  T  can  be  far  away  from  the  crack  tip  region  where  the  numerical 
evaluation  of  rj>,  wp,  dUj/dX^  and  oj’^1  is  usually  more  challenging. 
We  implemented  Eqs.  (30)  and  (31)  using  the  energy  domain 
integral  method  originally  proposed  by  Shih  et  al.  [40],  The  details 
of  this  numerical  implementation  can  be  found  in  Appendix  C. 

3.  Results  and  discussions 

The  analyses  carried  out  concern  the  energy  release  rate  for 
stationary  cracks.  Specifically,  a  large-deformation,  mixed  finite 
element  framework  of  the  rate  form  is  used  to  analyze  the  coupling 
between  diffusion  and  stress  development  [29,30],  Based  on 
numerical  solutions  obtained,  the  energy  release  rate  is  calculated 
using  the  path-independent  form  in  Section  2.  The  issues  of  focus 


1)  Flow  does  the  mechanical  driving  force  for  diffusion  affect 
crack-tip  fields  and  energy  release  rate? 

2)  Flow  does  plasticity  affect  the  energy  release  rate  in  the  pres¬ 
ence  of  diffusion/deformation  coupling?  and 

3)  How  does  lithiation-induced  softening  affect  fracture  and,  in 
turn,  battery  design? 

To  address  the  first  issue,  we  consider  a  highly-simplified  plane- 
strain  system  with  a  center  crack  in  a  2D  media  (Fig.  2a).  Lithiation- 
induced  softening  is  not  considered  in  this  case  and  the  material  is 
assumed  to  be  perfectly  elastic.  When  the  crack  is  loaded  by 
a  remotely  applied  stress  the  effect  of  diffusion/deformation 
coupling  can  be  analyzed,  especially  in  terms  of  crack-tip  fields  and 
energy  release  rate. 

Based  on  the  insights  gained  from  the  above  LEFM  problem, 
a  more  realistic  case  with  material  inelasticity  will  be  considered  to 
address  the  second  and  the  third  issues.  Specifically,  a  thin-film 
electrode  with  surface  cracks  during  delithiation  will  be  consid¬ 
ered  (cf.  Fig.  2b).  Analyses  in  this  regard  will  be  carried  out  in  two 
steps.  The  first  step  involves  constant  mechanical  properties  that 
are  independent  of  Li  concentration  (Section  3.2.1 ).  After  that,  the 
implications  of  lithiation-induced  softening  will  be  analyzed 
(Section  3.2.2). 

Depending  on  the  design  and  material,  stresses  in  an  electrode 
may  arise  via  three  mechanisms:  (1)  Li  concentration  inhomoge¬ 
neity  due  to  finite  diffusivity,  (2)  lattice  mismatch  at  phase 
boundaries,  and  (3)  constraining  by  external  agencies  such  as 
substrate  or  current  collector.  In  order  to  focus  on  the  effect  of 
deformation/diffusion  coupling  and  lithiation-induced  softening 
on  the  fracture  driving  forces,  we  only  consider  the  third  mecha¬ 
nism  (mechanical  constraint)  in  this  paper.  To  this  end,  we  assume 
that  the  initial  material  (Li/Si  in  this  paper)  is  fully  amorphized  and 
that  the  characteristic  loading  time  T0  is  much  longer  than  the 
characteristic  diffusion  time  tl,~H2/D ft  so  that  the  stresses  in  the 
electrode  are  entirely  due  to  external  loading  (for  Fig.  2a)  or 
mechanical  constraint  (for  Fig.  2b).  The  insight  thus  obtained  can  be 
easily  extended  to  situations  with  tl'~T0,  which  will  not  be  the 
focus  of  this  paper. 

3.1.  Linear  elastic  case 


Fig.  2.  (a)  A  highly-simplified  plane-strain  system  with  a  center  crack  in  a  2D  media, 
loaded  by  a  remote  stress.  This  configuration  will  be  used  to  analyze  the  effect  of 
diffusion/deformation  coupling  in  the  LEFM  regime,  (b)  A  thin-film  electrode  with 
periodic  surface  pre-cracks,  initial  homogeneous  Li  concentration  la,  and  zero  initial 
stress  undergoing  discharge  at  a  constant  surface  outflux.  The  pre-crack  is  modeled  as 
a  notch  with  a  small  but  finite  tip  radius  of  p  -c  a  for  numerical  stability. 


We  first  consider  the  linear  elastic  situation  such  that  the  yield 
stress  ffy  is  infinite  and  Young’s  modulus  E  and  Poisson  ratio  v  are 
2-independent.  The  geometry  of  the  stress-controlled  configura¬ 
tion  (cf.  Fig.  2a)  used  here  can  be  characterized  by  a  single 
parameter  a  which  measures  the  crack  length,  as  long  as  a.  The 


182 


Y.F.  Gao,  M.  Zhou  /  Journal  of  Power  Sources  230  (2013)  176-193 


Table  1 

Concentration-independent  material  parameters  for  Li^Si  used  in  the  LEFM  analysis. 


Parameter 

Symbol 

Typical  value 

Reference/ 

Partial  atomic  volume  of  Li 

qU(SF) 

14.2  A3 

[1] 

Si  concentration  in  fully-discharged 

co 

49.3  atoms  m 

m-3  [41] 

Young's  modulus 

£ 

80  GPa 

[42] 

Diffusivity 

Dq 

10-11  cm2  s“ 

1  [43-47] 

Poisson  ratio 

0.22 

[HI 

Maximum  charging  limit  of  LigSi 

fmax 

3.75 

[48] 

Chemical-to-mechanical  partial 

fl™ 

0.7 

fiUWc8 

Mechanical-to-chemical 
dimensionless  partial 
atomic  volume 

6mc 

274 

GU(SF)E/M 

pre-crack  is  modeled  as  a  notch  with  a  small  but  finite  tip  radius  of 
p  <c  a  for  numerical  stability.  Starting  from  a  stress-free  state  with 
homogeneous  initial  SOC  of  S0,  the  system  is  loaded  by  a  remote 
stress  of  rr0 „  =  (t/T0)E  with  xLi  -c  T0.  Under  the  condition  of 
zLl  -c  T0,  the  concentration  and  stress  fields,  hence  the  fracture 
driving  forces,  do  not  depend  on  Dg  explicitly.  Our  task  is  to  analyze 
the  effect  of  mechanical  driving  force  for  diffusion  on  crack-tip 
fields  and  the  energy  release  rate. 

Four  dimensionless  parameters  can  be  identified  here:  the 
charging  limit  ;max,  the  Poisson’s  ratio  v,  and  dimensionless  partial 
atomic  volumes  QCM  =  £?Ll(SF)Csi  and  QMC=Qu^SF)E/kB8.  The 
meaning  of  £2™  is  straightforward:  it  characterizes  the  volume 
expansion  ratio  when  Li  concentration  changes.  Specifically,  the 
ratio  between  the  volume  at  lithiated  state  LifH  and  the  volume  at 
fully-delithiated  state  Li0H  is  JSF  (?)  =  1  +  ?  £2™.  £3™,  therefore, 
determines  the  effect  of  chemical  diffusion  to  mechanical  defor¬ 
mation  (Chemical-to-Mechanical,  CM)  when  composition  changes. 
To  illustrate  the  significance  of  £2  ,  on  the  other  hand,  we  consider 

the  dimensionless  chemical  potential  pL,/kB8  which  can  be  easily 
shown  from  Eqs.  (A.13)-(A.15)  to  be 


kBe 


(32) 


The  first  term  of  Eq.  (32)  is  the  chemical  contribution  to  pLi/kB8, 
the  second  term  is  due  to  mechanical  stresses.  Since  om/E  is 
proportional  to  elastic  strain  ie  [cf.  Eq.  (A.l)  and  (A.2)],  one  can 
immediately  conclude  from  Fields  law  [Eq.  (A.10)]  that  £2  is  the 


parameter  which  controls  the  mechanical  driving  force  to  chemical 
diffusion  ( Mechanical- to-Chemi cal,  MC).  Specifically,  if  QMC  =  0 
there  would  be  no^stress  effect  on  J,J. 

The  role  of  £2  has  been  widely  studied,  especially  in  the 
context  of  thermally-induced  fracture.  The  significance  of  £2  and 
associated  effect  on  fracture,  however,  are  still  not  well- 
understood.  To  see  how  £2MC  affects  crack-tip  fields  and  energy 
release  rate,  we  used  the  typical  room-temperature  (0  =  300  K) 
values  of  material  parameters  for  LifSi  (H  -  Si  for  Li;H).  The 
...  ...  ~mc  amc 

discussions  compare  two  scenarios:  one  with  £2  —  £2LiSi  as  listed 

in  Table  1,  and  the  other  with  £2MC  =  0. 

3.1.1.  In-plane  stresses  and  stress  intensity  factor 

Fig.  3a  shows  the  redistribution  of  Li  concentration  S  =  ?/fmax 
when  an  initially  homogeneous  electrode  (with  a  =  500  nm  and 
S0  =  0.4)  is  loaded  by  a  remote  stress  of  <r«  =  0.005E,  under  the 
condition  of  £3  =  £3usi.  The  spatial  coordinates  are  normalized 

by  the  crack  length  a.  It  can  be  seen  that  2  deviates  from  the 
average  value  2  =  30  such  that  AE  =  E  -  20  >  0  in  front  of  the 
crack  tip.  This  redistribution,  driven  by  the  crack-tip  hydrostatic 
stress  <im,  in  turn  relaxes  om  (Fig.  3b),  as  shown  by  the  much  lower 
<rm  levels  for  the  fully-coupled  case  (£2  =  £2usi)  than  the  levels 

for  the  purely  mechanical  case  (£2MCm=  0). 

Given  the  dependence  of  am  on  £3  ,  one  may  conjecture  that  Li 

redistribution  lowers  all  components  of  crack-tip  stresses.  Quite 
counter-intuitively,  this  is  not  the  case.  As  shown  in  Fig.  4(a),  the 
in-plane  stress  components  and  ayy)  are  both  independent  of 

£3MC.  The  relaxation  of  am  is  entirely  due  to  the  change  in  <Tzz- 
-MC 

Indeed,  when  £2  =0  the  problem  reduces  to  an  elementary 

plane-strain  elasticity  problem  with  <rzz  =  v(<jxx  +  ayy)  >  0.  On  the 
other  hand  for  flMC  =  £2^,  a7Z<0  [cf.  Fig.  4(b)],  Such  a  tensile 
to  compressive  change  of  <rzz  significantly  increases  the  crack-tip 
von  Mises  stress  <rmises  [compare  Fig.  4(c)  and  (d)].  The  conse¬ 
quence  of  this  increase  of  amises  will  be  discussed  further  later 
(cf.  Fig.  8). 

In  terms  of  stress  intensity  factor  how  does  the  coupled  case 
here  differ  from  the  pure  linear  elastic  case  in  without  chemical 
transport  of  Li?  This  issue  warrants  a  careful  examination.  If  we 
follow  the  classic  linear  elastic  fracture  mechanics  (LEFM)  defini¬ 
tion,  the  single  parameter  I< i  that  determines  the  in-plane  stress 
fields  near  the  crack  tip  is  evaluated  through 


Fig.  3.  Distributions  of  Li  concentration  and  hydrostatic  stress  near  the  crack  tip  under  a  remote  stress  of  <r«,  =  0.005E.  (a)  deviation  of  concentration  from  its  spatial  average  when 
flMC  =  (b)  normalized  hydrostatic  stress  <rm/irx  along  they  =  0  line  for  the  fully-coupled  (fiMC  =  fl^si)  ar>d  the  purely  mechanical  (flMC  =  0)  cases.  The  inset  in  (b)  shows 

the  corresponding  om/ax  contours.  The  length  of  the  center-crack  is  a  =  500  nm  and  the  initial  SOC  is  E0  =  S  =  0.4.  All  spatial  dimensions  are  normalized  by  crack  length  a. 
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Fig.  4.  Distributions  of  ir**,  <ryy,  oiz  and  the  von  Mises  effective  stress  ffmjsCS  normalized  by  the  remote  stress  <r„  =  0.005E.  (a-b)  Normalized  in-plane  stresses  (a)  and  out-of-plane 
stress  (b)  along  y  =  0.  (c-d)  Normalized  von  Mises  stress  for  fully-coupled  (c)  and  purely  mechanical  (d)  cases.  The  length  of  the  center-crack  is  500  nm  and  the  initial  SOC  is 
S0  =  2  =  0.4.  All  the  spatial  dimensions  are  normalized  by  the  crack  length  a. 


Er»  -  JgjlM-  «J  -  t.Z).  (33) 

where  (r,  6)  are  the  polar  coordinates  with  origin  at  the  crack  tip. 
For  both  plane-strain  and  plane-stress  conditions,  Eq.  (33)  has  the 
same  form  of  fj{d).  Effect  of  out-of-plane  constraint  kicks  in  only 
through  out-of-plane  stress:  for  plane-strain,  0&  =  v(tjxx  +  <%); 
for  plane-stress  azz  =  0. 

Although  the  stress  states  in  Fig.  4  are  computed  under 
plane-strain  conditions,  the  relationship  of  (Jzz  =  v(dxx  +  fyy)  is 
no-longer  valid.  Instead,  ffzz  <  0  for  QMC  =  as  seen  in  Fig.  4b. 
Still,  the  results  in  Fig.  4a  indicate  that  the  in-plane  stresses  are 
independent  of  QMC.  In  other  words,  as  long  as  the  material  is 
elastic  and  sufficient  time  is  given  to  allow  diffusion  to  occur  (i.e. 
r0W),  the  in-plane  stresses  can  always  be  characterized  by  Eq. 
(33)  even  though  stresses  affect  Li  redistribution.  The  outcome  is 


that  K[  is  independent  of  the  accumulation  of  Li  at  the  crack  tip. 
Consequently,  Kj  is  incapable  of  capturing  the  full  effect  of 
diffusion-deformation  coupling  on  fracture. 

3.1.2.  Energy  release  rate 

Similar  findings  on  ffi  have  been  made  for  hydrogen  embrittle¬ 
ment  of  metals  [27]  and  transformation  toughening  of  ceramic 
materials  [49],  in  the  sense  that  the  in-plane  stresses  are  indepen¬ 
dent  of  dilatational  eigen-strains  near  the  stationary  crack  tip.  This 
conclusion  on  in-plane  stresses  can  also  be  proved  analytically  using 
the  linear  elastic  compatibility  conditions  (which  will  not  be 
elaborated  here).  However,  previous  studies  have  not  clarified  if  the 
energy  release  rate,  a  more  universal  parameter  for  fracture  analysis, 
changes  as  Li  redistributes.  We  provide  an  analysis  in  this  regard. 

Fig.  5a  shows  the  calculated  J  values  for  the  range  of  initial 
concentration  of  0.2  <  30  <  0.999.  The  circular  symbols  denote  the 
non-coupled  energy  release  rate  yNCP=J(f2MC  =  0)  and  the  solid 


a„/E  x  10'3  s0 

Fig.  5.  Path-independent  J  integral  for  a  stationary  crack  in  Fig.  2a  using  Eq.  (31 ).  (a)  Energy  release  rate  for  the  fully-coupled  case ]CP=J(QMC  =  -Q^  jas  a  function  of  remote  stress 
for  different  levels  of  initial  concentration  20,  the  circular  symbols  denote  JNCP=J(Q  =  0)  as  the  baseline  for  comparison,  (b)  Dependence  of  f-p  on  initial  concentration  when  Kt 
(thus/®3”)  is  fixed  at  a  value  that  corresponds  to  a  =  500  nm  and  =  0.005E. 
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the  crack  tip  is  subject  to  uniaxial  tensile  stresses,  Li  transported  from  a  reservoir  may 
overcome  the  elastic  Poisson’s  effect  and  induce  an  overall  expansion  in  the  transverse 
directions.  The  cuboid  with  dotted  boundary  denotes  the  initial  configuration  of  the 
isolated  near-tip  system  before  stretching,  and  the  cuboid  with  solid  boundary 
represent  its  stretched  counterpart. 

lines  represent  the  fully-coupled  values  ofJcp=J(QMC  =  QyS^).  For 
QMC  =  0JNCP  coincides  with  the  purely  mechanical  LEFM  result  of 
JNCP  /(aE)  —  7t(l  -  v2)(aoo/E)2.  For  QMC  —  two  observations 
can  be  made.  First,  JCP  depends  on  S0;  second,  JCP  >  JNCP. 
Specifically,  for  all  initial  concentrations  Jcp  is  always  higher  than 
JNCP,  except  for  B0  —  0.999  when  the  two  are  essentially  equal.  This 
trend  is  more  clearly  seen  in  Fig.  5b,  which  shows  JCP//NCP  as 
a  function  of  So  (crack  length  of  500  nm  and  remote  load 
ff„„  =  0.005F). 

The  S0 -dependence  of  Jcp  is  due  to  the  fact  that  Li  content  S 
cannot  exceed  1,  as  embodied  by  the  fact  that  nu/kB0->  <*>  when 
S— >1  [cf.  Eq.  (32)].  If  the  initial  content  S0  is  close  to  1,  redistri¬ 
bution  of  Li  cannot  happen  and  Jcp  approaches  JNCP.  The  plateau  for 
S0  <  0.8  [cf.  Fig.  3b],  on  the  other  hand,  is  due  to  the  fact  that  the 
relaxation  of  <jm  is  nearly  complete  such  that  changes  in  S0  in  the 
range  of  0.2  <  S0  <  0.8  bear  no  significant  influence  on  Jcp. 

The  finding  that  JCP/JNCP  >  1  appears  counterintuitive  at  the  first 
glance:  since  Li  redistribution  is  a  stress  relaxation  mechanism, 
why  would  this  redistribution  cause  the  energy  release  rate  to 
increase?  It  turns  out  that  this  anomaly  can  be  associated  with 
a  negative  effective  Poisson’s  ratio  effect  near  the  crack  tip.  If  we 
conceptually  isolate  the  LifSi  material  near  the  crack  tip,  the  bulk  of 
LifSi  far  from  the  crack  can  be  regarded  as  a  reservoir  with  fixed  Li 
chemical  potential  fi^s  [Fig.  6],  Consider  a  uniaxial  stretch  applied 
to  the  near-tip  material.  During  the  stretch,  fi%p  decreases  as  the 
hydrostatic  stress  becomes  tensile.  The  associated  chemical 
potential  difference  fi^p  -  fi^es  <  0  drives  Li  from  the  reservoir  into 
the  crack  tip.  If  the  amount  of  Li  so  transported  is  large  enough,  the 
associated  lithiation-induced  volume  change  may  overcome  the 
elastic  Poisson  effect  and  lead  to  transverse  expansion.  The 
outcome  is  effectively  a  negative  Poisson’s  ratio  —  tensile  loading 
causes  lateral  dimensions  to  increase.  This  is  exactly  what  happens 
here  for  LifSi. 

To  quantify  this  effect,  we  consider  the  limiting  case  with 
Djj  =  oo,  Q  =oo  and  ?max  -  Under  this  condition,  the 
chemical  potential  in  the  entire  material  is  always  in  equilibrium 
and  the  mechanical  contribution  to  diffusion  -QMC(<rm/E)  domi¬ 
nates  the  chemical  contribution  $ / (Usd)  +  ln[S/(l  -S)]  [cf.  Eq. 
(32)].  Consequently,  om  in  the  near-tip  region  is  pinned  at  the  same 
level  as  that  in  the  remote  regions. 

-MC 

For  Q  =  0,  the  elastic  response  of  the  material  is  character¬ 
ized  by  bulk  modulus  K  and  shear  modulus  G.  For  Q  =  oo ,  on  the 
other  hand,  the  overall  deformation  is  governed  by  effective  moduli 
Geff  =  G,  Keff  —  0,  and  veS  =  -1,  reflecting  the  fact  that  near-tip 
om  is  pined  at  the  remote  level. 


Now  that  k|  is  independent  of  QUC  [Discussions  in  3.1.1  ],  the 
introduction  of  Ge %  and  KeS  allows  J(QMC  =  oo)  to  be  calculated 
from  K i  using  the  classic  LEFM  relation.  Specifically, 


j{pUC 


1  -  (irfff)2 
£<# 


Kj 


4 

Geff  G  ' 


(34) 


In  contrast,  the  energy  release  rate  for  the  purely  mechanical 
case  is 


JNCP  =  j(qMC  =  0)  =l^Kf=l^Kf.  (35) 

We  therefore  have  shown  that 


K**--)  2 

j(£mc  =  o)  i- 


(36) 


Indeed,  for  v  =  0.22  the  ratio  given  by  Eq.  (36)  is  2.56,  which  is 
fairly  close  to  the  plateau  value  (  =  2.2)  of  JCP/JNCP  in  Fig.  5b.  The 
difference  is  due  to  the  fact  that  the  ideal  case  of  i2MC  =  °°  [Eq. 
(36)]  overestimates  the  relaxation  of  om  compared  with  the 
realistic  case  of  fiMC  =  [Fig.  5b], 

From  an  energetic  point  of  view,  Jcp  >  JNCP  embodies  the  fact 
that  the  redistribution  of  Li  provides  another  source  of  energy  for 
the  fracture  driving  force  besides  the  mechanical  fields.  Two 
complementary  processes  make  this  additional  source  of  energy 
available  for  fracture.  The  first  is  energy  storage.  For  QMCj=  0,  the 
external  agency  (<r „  here)  must  provide  extra  work  compared  with 
the  case  with  Q  -  0  in  order  to  induce  the  same  remote  fields. 
This  extra  work  is  stored  as  latent  energy  in  the  form  of  inhomo¬ 
geneous  Li  concentration.  The  second  mechanism  is  energy 
retrieval.  Through  the  agency  of  composition-change-induced 
stresses,  the  latent  energy  is  converted  back  into  mechanical 
energy,  thereby  increasing  the  fracture  driving  force. 

The  storage  process  is  controlled  by  Q  and  the  retrieval 

ii  j  i  ;aCM  .  ~  MC  -CM  . 

process  is  controlled  by  Q  .If  either  Q  or  Q  is  zero,  the 
coupling  is  severed  and  the  energy  release  rate  reduces  to  the 
purely  mechanical  level.  On  the  other  hand,  the  analysis  in  Eqs. 
(34)— (36)  provides  an  upper  bound  for  the  change  of  energy 
release  rate  due  to  this  storage— retrieval  mechanism.  This  limit  is 
reached  when  the  hydrostatic  stress  field  is  completely  relaxed  (to 
the  remote  level)  so  that  full  redistribution  of  Li  has  occurred.  This 
upper  bound  is  universal  or  the  same  for  all  cases  regardless  of 
material  properties  and  loading,  as  long  as  the  material  response  is 
elastic  and  the  loading  is  slow. 

It  should  be  noted  that  the  anti-shielding  effect  of  J  due  to 
deformation-diffusion  coupling  (i.e.,  the  fact  that  Jcp  >  JNCP)  in  this 
paper  is  evaluated  for  the  steady-state  of  Li  distribution  such  that 
T0»Til.  Under  this  condition,  stresses  arise  solely  because  of 
external  constraint  [Fig.  2b]  or  mechanical  loading  [Fig.  2a],  As 
pointed  out  earlier,  in  general  stresses  may  also  develop  during 
charge/discharge  due  to  concentration  inhomogeneity  when 
tu~T0,  as  seen  in  free-standing  amorphous  LiSi  particles  or  NWs 
[16,19,50,51],  For  such  situations  with  transient  effect, 
deformation-diffusion  full  coupling  would  manifest  itself  in  two 
aspects.  First,  accumulation  of  Li  near  the  crack  tip,  again,  leads  to 
an  anti-shielding  effect  on  J  which  increases  the  magnitude  of  the 
energy  release  rate.  Second,  remote  concentrations  and  stresses  far 
from  the  crack  are  governed  by  an  effective  diffusivity  Deg  [11,16], 
This  Deff  leads  to  a  stress-enhanced-diffusion  (SED)  effect  that 
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lowers  the  remote  stresses  [16]  and  hence  J.  In  other  words,  the  first 
effect  increases  J  while  the  second  effect  decreases  J.  The  final 
outcome  is  the  combined  effect  of  the  two.  Depending  on  the  ratio 
of  a/H,  the  characteristic  time  scales  for  the  two  effects  can  be 
different.  This  transient  situation  of  rL'  ~  T0  is  not  the  focus  of  this 
paper.  Instead,  configurations  shown  in  Fig.  2  with  T0^>tu  are 
considered  in  order  to  isolate  the  effect  of  Li  accumulation  near  the 
crack  tip  (the  first  effect).  The  insight  thus  obtained  can  be  easily 
combined  with  the  effect  of  SED  for  far  fields  and  thereby  extended 
to  situations  with  tu~T0. 

3.1.3.  Remarks  on  IQ-based  and  J -based  fracture  approaches 

We  have  shown  earlier  that  IQ  is  independent  of  Q  ,  i.e.  Ki  is 
unaffected  by  the  deformation/diffusion  coupling  in  the  long-term 
limit  (i.e.  T0 s> zu).  In  contrast,  Fig.  5  and  Eqs.  (34)— (36)  show 
clearly  that  the  coupled  energy  release  rate  does  depend  on  QMC . 
Even  when  the  conditions  of  plane  strain  are  maintained,  the 
relationship  between  J  and  IQ  is  no  longer  unique,  in  the  sense  that 
the  coupled  energy  release  depends  on  Q  ,  Q  and  E0.  Given  this 
dichotomy,  we  ponder  the  question  of  which  criterion,  K\  —  IQC  or 
J  =  Jcr .  is  more  appropriate. 

To  answer  this  question,  we  first  note  that  a  similar  dichotomy 
exists  in  classic  fracture  mechanics  in  that  the  relationship  between 
J  and  Ki  is  not  always  unique  there  either.  Specifically,  under  the 
same  in-plane  loading,  JNCP  =  (1  -  v2)Kf/E  for  plane-strain  and 
jNCP  _  i(2 /p  for  plane-stress.  This  difference  arises  due  to  the 
difference  in  the  out-of-plane  constraint  and,  hence,  crack  tip 
stress  triaxiality.  Here,  however,  the  situation  is  more 
complicated  because  the  crack-tip  stress  triaxiality  is  not  only 
determined  by  out-of-plane  constraint  bjrj^alsgJjy  the  redistribu¬ 
tion  of  Li.  Depending  on  the  values  of  Q  ,  Q  and  S0,  azz  may 
assume  any  value  between  r(<r»H- <tyy)  and  -{<Jxx  +  ayy),  and  J 
may  be  of  any  value  between  (1  —v2)Kf/E  [for  QMC  —  0]  and 
Kf/G  =  2(1  +  v)Kf/E  [for  fiMC  =  »],  even  under  the  same  plane- 
strain  conditions. 

In  the  pure  mechanical  case,  the  fracture  criterion  can  be 
stated  as 

*(*"*)  =  ^critical  (dcrac*£) ,  Or 

j(dcrack)  =  Jcr  (dcrack) . 

The  effect  of  crack-tip  load  triaxiality  can  be  lumped  into  the 
thickness-dependence  of  fracture  toughness  KCriticai(dcrack)  which  is 
higher  than  the  plane  strain  fracture  toughness  K\c  and  lower  than 
the  plane  stress  fracture  toughness  IQ  such  that  Klc  <  Kcriticai(dcrack) 
<  Kc.  This  handling  allows  the  stress  intensity  factor  IQ  calculated 
solely  from  in-plane  stresses,  which  do  not  depend  on  out-of-plane 
constraints,  to  be  compared  with  the  thickness-dependent  fracture 
toughness  Kcritical(dcrack)  to  determine  the  onset  of  fracture.  As 
a  conservative  measure,  the  plane-strain  fracture  toughness  IQC  is 
most  often  used.  A  similar  treatment  is  reflected  in  the  ./-integral 
based  fracture  criterion  above  in  that  both  J  and  Jcr  are  thickness- 
dependent. 

Similarly,  it  is  also  possible  to  state  the  fracture  criterion  for 
electrode  materials  which  deforms  only  elastically  and  experience 
diffusion/deformation  coupling  as 

K(cF^)  =  Kcriticai (dcradi,  20) ,  or  1 
j(dcrack,  S0)  =  JCR  (dcrack,  20) ,  j 


as  long  as  both  sides  of  these  relations  are  specific  to  the  same 
thickness  and,  in  addition,  the  right  hand  sides  are  also  specific  to 
the  initial  state  of  charge  B0.  At  this  stage,  it  is  unclear  under  what 
conditions  the  critical  value  functions  Kcriticai(dcrack,E 0)  and 
JcR(dcrack ,  B0)  would  assume  their  minima  or  most  conservative 
values.  These  failure  envelopes  separating  the  safe  and  unsafe 
regions  in  the  (dcrack,  S0)  space  can  only  be  determined  by 
experiments. 

As  a  concluding  remark  to  this  section,  it  should  be  noted  that 
the  results  in  Figs.  (3)— (5)  are  based  on  the  highly-idealized 
assumption  that  the  material  is  linear  elastic.  Understandably,  the 
behavior  of  real  Li/Si,  which  is  capable  of  deforming  plastically,  is 
different.  The  results  in  this  section  may  be  more  applicable  to 
cathode  materials  such  as  LiFePC>4  which  is  generally  quite  brittle. 
When  significant  (large-scale)  plasticity  is  involved,  the  K-based 
criterion  [the  first  of  Eq.  (38)]  should  not  be  used  and  the  ./-based 
criterion  [the  second  of  Eq.  (38)]  is  the  only  sensible  approach. 
Nevertheless,  the  analysis  in  this  section  yields  insight  into  the 
behavior  of  alloy-based  electrode  materials.  This  discussion  lays  the 
foundation  for  further  discussions  accounting  the  effect  of  plasticity 
and  lithiation-induced  softening  in  the  following  sections. 

3.2.  Elasto-plastic  case 

We  now  turn  our  attention  to  the  situation  when  the  electrode 
material  deforms  plastically,  especially  when  plastic  deformation  is 
large  enough  to  render  the  LEFM  approach  using  K  invalid.  To  this 
end,  we  first  consider  a  simplified  case  such  that  the  yield  stress  aY 
is  assumed  to  be  independent  of  B0  (Section  3.2.1 ).  For  this  purpose, 
we  assume  the  same  material  properties  and  dimensionless 
parameters  as  listed  in  Table  1,  with  only  one  extra  dimensionless 
number  #y=%/E  <  1  which  quantifies  concentration- 
independent  yield  stress  ay.  The  case  with  concentration- 
dependent  aY  and  elastic  moduli  will  be  considered  in  Section  3.2.2. 

When  the  material  deforms  plastically,  the  stress-controlled 
configuration  [Fig.  2(a)]  used  to  facilitate  comparison  with  pure 
mechanical  LEFM  solutions  is  no  longer  appropriate  because  any 
remote  stress  a „  exceeding  aY  would  lead  to  large-scale  plastic 
deformation.  In  this  section,  we  therefore  use  the  configuration 
shown  in  Fig.  2b  which  involves  a  thin-film  electrode  with  initial 
thickness  H  and  SOC  S0  under  galvanostatic  discharging.  Again,  the 
pre-defect  is  a  surface  crack  with  length  a  and  a  small  but  finite 
crack-tip  radius  of  p  <  a.  The  initial  stress-free  state  (with  SOC  of 
S0)  is  taken  to  be  the  Lagrangian  configuration  in  which  the  energy 
release  rate  is  calculated  via  Eq.  (31).  A  surface  outflux  of  Jsurf  is 
specified  to  simulate  galvanostatic  discharging  such  that  the 
average  SOC  S  in  the  film  decreases  with  time  t  according  to  S(t)  = 
S0  -  t/To  for  t  <  Eq/Tq,  where  T0  =  (1  hour)/(C  -  rate)  is  the 
nominal  discharge  time.  The  nominal  thickness  H0  is  the  thickness 
of  the  film  in  the  fully  discharged  (i.e.,  3  =  0  throughout  the 
electrode)  state.  It  is  related  to  the  initial  thickness  H  via 
H~Jsf(E0)H0.  The  approximation  is  associated  with  elastic  strains 
and  therefore  is  very  reasonable.  Unless  indicated  otherwise,  H0  is 
used  because  it  is  more  convenient  and  almost  exclusively  used  in 
the  literature  [20,52—54], 

For  the  calculation  in  this  section,  a  typical  film  thickness  of 
Hq  =  200  nm  and  nominal  discharging  time  of  T0  =  50  h  are 
considered.  These  values  are  such  that  T0  » ih  and  are  in  the  range 
of  available  experimental  data.  In  addition,  a/H  =  a/\JSF(E0)H0]  < 
1  /4,  so  that  the  film/substrate  interface  does  not  significantly  affect 
crack-tip  fields. 

3.2.1.  Effect  of  global  yielding  on  energy  release  rate 

For  galvanostatic  discharging,  the  dimensionless  discharge  level 
AE=E0  -  E(t)  =  t/To  can  be  used  to  conveniently  characterize  the 


(37) 
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amount  of  Li  extracted  during  the  process.  When  T0 ;» iLl,  all 
stresses  are  due  to  constraint  at  the  film-substrate  interface,  and 
the  bi-axial  stress  state  far  from  the  crack  can  be  written  as 


„  .jt/W 

<7 XX  =  (Tyy  =  <Ty, 


for  AS  <  A3y, 
for  AS  >  ASy, 


(39) 


where  ASy  is  the  discharge  level  at  which  axx  and  ayy  reach  ay  so 
that  the  film  begins  to  yield  globally.  ASy  is  approximately  inde¬ 
pendent  of  QMC  so  long  as  T0  »  xLl,  because  flMC  only  affects  the  Li 
distribution  near  the  crack  and  has  no  bearings  on  stresses  far  away 
from  the  crack  (y»H)  where  S  —  E0  -  ASy  is  homogeneous. 

Fig.  7  shows  the  distribution  of  effective  plastic  strain 
ipeff=  Jo  y(2/3)DP  :  DPdt  for  a  film  with  aY  =  0.01  ,H0  =  200  nm, 
a  =  50  nm  and  S0  =  0.5,  at  AS  =  ASy  as  specified  by  Eq.  (39).  The 
plastic  zone  indeed  extends  through  the  thickness  (horizontal)  of 
the  film  for  both  QMC  —  and  QMC  —  0,  confirming  the  onset  of 
full-scale  yielding.  The  remote  regions  show  slightly  later  yielding 
for  fiMC  =  flyjj  than  for  QMC  —  0,  mainly  because  of  the  finite 
dimension  in  the  y-direction.  This  slight  difference  does  not  affect 
the  main  conclusions  to  be  drawn.  On  the  other  hand,  the  plastic 

~  mc  ~  MC  ~  MC 

strain  near  the  crack  is  higher  for  Q  =  £?Lisi  than  for  Q  -  0,  in 
consistency  with  the  previous  finding  that  Li  redistribution 
increases  the  crack-tip  von  Mises  stress  [cf.  Fig.  4c  and  d].  This 
promotion  of  plastic  deformation,  together  with  the  redistribution 
of  Li,  is  also  responsible  for  the  observation  that  the  crack  opening  is 
wider  for  QMC  —  fiys,  than  for  fiMC  =  0  [see  Fig.  8a  and  b]. 

The  redistribution  of  Li  and  resultant  relaxation  of  am  at  ASy  are 
shown  in  Fig.  8.  Although  significant  plastic  strains  are  involved,  the 
redistribution  of  Li  [Fig.  8c]  is  rather  similar  to  that  in  the  fully 
elastic  case  [Fig.  3a[.  Again,  this  redistribution  causes  am  to  relax 
[Fig.  8a  and  b].  In  contrast  to  what  is  seen  previously,  the  in-plane 
stresses  are  no  longer  independent  of  Q  .  This  difference  indicates 
that  J  is  the  only  sensible  parameter  for  characterizing  the  fracture 
behavior.  It  turns  out  that  the  path-independent  form  of  J  in  Eq.  (31 ) 


Fig.  7.  Effective  plastic  strain  f,,  for  a  st< 
nominal  thickness  H0  =  200  nm  and  ii 
discharge  reaches  a  level  ASY  such  that  the  film  be 
cry  =  0.01.  (a)  for  fully-coupled  case  fiMC  =  (b)  £?jwitl1  ^MC  =  °- 


Fig.  8.  Distributions  of  Li  concentration  and  hydrostatic  stress  for  a  stationary  crack 
with  a  =  50  nm  in  a  film  with  nominal  thickness  H0  =  200  nm  and  initial 
concentration  S0  =  0.5.  The  condition  shown  corresponds  to  onset  of  global  yielding 
at  AS  =  ASy  with  crY  =  0.01.  (a)  hydrostatic  stress  normalized  by  yield  stress  for 
flMC  =  flysj,  (b)  hydrostatic  stress  normalized  by  yield  stress  for  fiMC  =  0,  and  (c) 
deviation  of  Li  concentration  from  its  average  value  3  =  30  -  ASy  for  flMC  = 

is  highly  robust  numerically,  with  less  than  0.5%  difference  between 
J  calculated  from  difference  contours  even  when  AS  reaches  ASy.  A 
brief  discussion  on  the  numerical  stability  and  errors  associated 
with  this  method  can  be  found  in  Appendix  D. 

Fig.  9  shows  the  evolution  of  J  [normalized  by  (Fa)]  for 
ay  =  0.01  (solid  black  lines).  Both  Ja‘  —  J(flMC  =  Qy^)  and JNCP  = 
J{Q  =  0)  are  plotted,  with  the  corresponding  curves  for  the 
purely  elastic  cases  (red  dashed  lines)  given  for  comparison.  Clearly, 
plasticity  significantly  affects  the  energy  release  rate  (both  Jcp  and 
JNCP),  except  for  the  very  initial  stages  of  discharge  (AS  <  0.01) 
when  the  plastic  zone  is  small  enough  so  that  LEFM  provides  a  good 
approximation. 
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Fig.  9.  Evolution  of  energy  release  rate  J  for  a  stationary  surface  crack  with  a  =  50  nm 
in  a  film  with  nominal  thickness  H0  =  200  nm  and  initial  concentration  S0  =  0.5. 
Solid  black  lines  indicate  normalized  values  of  Jcp  =  ](QMC  =  and 


Fig.  10.  Dependence  of  energy  release  on  yield  stress  aY.  The  length  of  the  surface 
crack  is  a  =  50  nm  and  the  nominal  film  thickness  is  H0  =  200  nm.  The  initial 
concentration  is  S0  =  0.5.  Full  deformation-diffusion  coupling  ( Q  =  flusl)  is 
considered. 


jncp  _  j(QMC  =  0)  for  <iy  =  0.01.  Dashed  red  lines  show  the  corresponding  Jcp  and 
f*  values  for  the  corresponding  elastic  cases.  (For  interpretation  of  the  references  to 
colour  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 

As  AE  increases,  both  Jcp  and  JNCP  begin  to  deviate  from  their 
LEFM  values.  The  most  striking  feature  is  the  sudden  slope  changes 
in  the  J  -  AE  curves  (for  both  fully-coupled  case  of  Q  —  Qusi  and 
purely  mechanical  case  of  QMC  —  0)  when  AS  reaches  the  global 


The  discharge  level  AECR  when  critical  condition  J  -  JCR  is 
reached  determines  the  maximum  extractable  amount  of  Li 
without  electrode  failure.  Designs  and  cycling  regimens  that 
maximize  AEcr  would  increase  the  utilizable  capacity  of  the  elec¬ 
trode  without  jeopardizing  its  cyclability.  For  the  same  JCR  and  a. 
Fig.  10  shows  that  choosing  materials  with  smaller  E  and  cry  (hence 
cry)  is  an  effective  way  to  avoid  fracture. 


yielding  threshold  AEY  [cf.  Fig.  7],  Before  AEY,  plastic  zone  size  ry  is 
small  relative  to  the  system  dimension,  and  the  effective  K i  and  J  can 
be  estimated  using  Irwin’s  effective  crack  length  aeg  =  a  +  ry  [55], 
Beyond  AE  =  AEY,  however,  the  film  yields  globally  and  Irwin’s 
assumption  no  longer  applies.  Instead,  plasticity  in  regions  far  away 
from  the  crack  tip  becomes  important.  This  prevents  the  remote 
stresses  from  increasing  further  and  leads  to  much  flatter  JCF  and 
JNCP  curves.  Such  global  yielding  is  important  because  it  signifi¬ 
cantly  reduces  the  energy  release  rates,  thereby  lowering  the 
fracture  tendency. 

The  shapes  of  the  J  -  AE  curves  in  Fig.  9  are  consistent  with  that 
was  found  by  Aoki  et  al.  [56]  who  calculated  the  energy  release  rate 
for  a  stationary  crack  in  a  thermally-loaded  plate.  What  they 
studied  corresponds  to  the  special  situation  with  Q  >0  and 
Q  =  0.  Our  calculation  confirms  that  the  overall  characteristics 
of  the  J  -  AS  curve,  especially  the  existence  of  inflection  points, 
-MC 

remain  even  for  fully-coupled  case  with  Q  =7=0.  More  importantly, 
Fig.  9  indicates  that  Jcl7jNCP  >  1  even  when  the  combined  effect  of 
plasticity  and  diffusion/deformation  coupling  is  considered, 
although  the  upper  bound  given  by  the  linear  elastic  solution  [Eqs. 
(34)-(36)]  can  no  longer  be  used  directly. 

Fig.  10  shows  the  dependence  of  J  on  rrY  for  the  fully-coupled 
case  of  QMC  —  For  both  <rY  =  0.005  and  cry  =  0.01,  an 
inflection  point  is  seen  in  the  J-AE  relation.  The  location  of  the 
inflection  point  for  each  case  corresponds  to  its  global-yielding 
threshold  AEY.  Overall,  the  reduction  in  ay  lowers  the  fracture 
tendency  by  lowering/ [e.g.J(<7y  =  0.005)  <  J(aY  —  0.01)  <  Jeiastic 
for  AE  >  0.05],  except  for  early  stages  before  global  yielding  occurs 
so  that  the  development  of  J  is  approximately  governed  by  the 
effective  crack  length  aejj  >  a.  In  terms  of  elastic  modulus,  note  that 
J  is  normalized  by  Ea  so  a  reduction  in  E  would  also  lead  to 
a  reduction  in  the  fracture  driving  forces. 


3.2.2.  Effect  of  lithiation-induced  softening 

What  we  learned  from  Fig.  10  can  be  important  for  alloy-based 
electrode  materials,  because  the  elastic  moduli  and  yield  stress  of 
such  materials  show  significant  dependence  on  Li  concentration.  If 
we  manage  to  find  a  composition  window  Ee  (E0  -  AECR,  E0)  for 
which  the  combination  of  elastic  modulus  and  yield  stresses 
minimizes  J,  we  can  operate  the  battery  in  such  a  window  to 
maximize  the  operational  capacity  AECR. 

To  find  such  a  window  for  Li/Si,  we  now  put  all  the 
insights  learned  thus  far  together  and  assess  the  combined 
effect  of  diffusion/deformation  coupling,  inelasticity,  and 
lithiation-induced  softening  on  fracturing  driving  force.  To  this 
end,  calculations  involving  composition-dependent  bulk 
modulus  K  =  (12.46?  +  65.44)/(l  +  ?)  GPa,  shear  modulus 
G  =  (7.63?  +  35.51)/(l+?)  GPa  [23],  and  yield  stress 
<jy  =  (1.75  +  0.167?)/(1  +  ?)  GPa,  are  conducted.  The  variation  of 
yield  stress  with  concentration  represents  a  simple  interpolation 
between  aY  =  1.75  GPa  at  ?  =  0  and  <rY  =  0.5  GPa  at  ?  —  fmax 
[57,58]  using  the  rule  of  mixture.  Although  it  is  likely  that  JCR  for  Li/ 
Si  is  composition-dependent,  we  adopt  a  constant  value  of 
JcR  -  3  J/m2  due  to  lack  of  experimental  data.  This  critical  value  is 
estimated  using  the  toughness  for  pure  silicon  which  is  typically 
between  2.5~4.3  J/m2  [59]. 

Since  the  mechanical  properties  here  are  concentration- 
dependent,  dimensionless  parameters  QMC  and  dY  are  no  longer 
convenient,  and  all  quantities  will  be  stated  in  their  respective 
physical  units.  The  system  configuration  and  discharge  rate  are  the 
same  as  those  in  Section  3.2.1  [Fig.  2b],  which  is  relevant  for  thin- 
film  electrodes  under  galvanostatic  discharge. 

Fig.  11  shows  the  dependence  of  J  development  on  initial 
concentration  E0,  for  a  film  with  nominal  thickness  of 
H0  =  200  nm  and  pre-crack  length  of  a  =  20  nm.  For  each  J  -  AE 
curve,  the  critical  discharge  level  AEcr  at  which  J  —  JCR  is  marked. 
Two  scenarios  can  be  identified.  For  lower  E0  [e.g.  S0  <  0.6]  the 
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Fig.  11.  Dependence  of  energy  release  rate  on  discharge  level  AS  and  initial  concen¬ 
tration  S0  for  a  crack  length  of  20  nm.  Full  deformation-diffusion  coupling  is 


energy  release  rate  reaches  its  critical  value  before  the  onset  of 
full-scale  yielding.  For  higher  30  [e.g.  S0  >  0.8],  on  the  other  hand, 
the  material  is  softer  and  the  inflection  point  AEy  is  reached  before 
fracture. 

In  the  first  scenario,  fracture  is  essentially  governed  by  Irwin’s 
effective  crack  length.  Under  this  condition,  the  A 3CR  -  S0 
dependency  is  mainly  due  to  lithiation-induced  elastic  softening 
(i.e.  reduction  of  elastic  moduli  as  Li  concentration  increases). 
Specifically,  since  J  is  approximately  proportional  to  the  Young’s 
modulus  [note  that  J  is  normalized  by  ( Ea )  in  Fig.  10],  the  fracture 
driving  force  is  lower  when  the  operation  window  is  chosen  at 
higher  concentrations.  In  the  second  scenario,  relaxation  due  to 
global  yielding  comes  into  play,  leading  to  much  wider  windows  of 
Eg  (E0  -  A Equ,  E0 )  than  that  in  the  first  scenario.  The  increase  in 
AECr  due  to  global  yielding  is  much  more  significant  compared 
with  that  due  to  elastic  softening.  For  all  cases,  the  maximum 
allowed  operational  capacity  AECr  is  always  higher  when  30  is 
higher,  as  a  combined  consequence  of  the  lithiation-induced 
softening  of  E  and  oy. 

We  conclude  our  discussion  with  a  design  map  [Fig.  12] 
which  shows  the  maximum  utilizable  amount  of  Li  (as  measured  by 


Fig.  12.  Design  map  showing  the  maximum  extractable  amount  of  lithium  AS  without 
crack  growth  as  a  function  of  initial  SOC  S0  for  different  pre-crack  sizes.  Full  defor¬ 
mation-diffusion  coupling  is  considered. 


A3)  as  a  function  of  S0  when  pre-cracks  with  different  sizes  are 
present.  Again,  Eg  (30  -  AECr,  S0)  is  the  safe  window  in  which  the 
electrode  can  be  operated  without  fracture.  For  each  pre-crack 
length  a,  the  corresponding  AECR  -  30  curve  delineates  the 
boundary  between  the  unsafe  region  (upper  left)  and  the  safe 
region  (bottom  right).  In  real  thin-film  electrodes,  the  smoothness 
of  electrode  surfaces  is  always  limited  by  the  manufacturing 
and  operating  conditions,  and  the  existence  of  surface  defects  is 
inevitable.  Fig.  12  indicates  that  for  the  same  operational 
capacity  (AS),  a  thin-film  electrode  is  always  more  defect-tolerant 
when  the  changing/discharging  regimen  is  designed  such  that 
the  battery  operates  at  higher  concentration  windows  of 
Ee(E0  -  AEcr,  E0).  _ 

As  a  final  remark,  it  should  be  noted  that  AEcft  in  Figs.  11  and  12 
is  calculated  based  on  the  assumption  of  a  constant  JCr.  For  real 
alloy-based  electrode  materials,  it  is  very  likely  that  JCr  is 
concentration-dependent.  Taking  conventional  engineering  mate¬ 
rial  as  reference,  the  general  trend  for  fracture  toughness  is  that  Kic 
is  higher  for  materials  with  lower  try.  If  this  trend  is  also  true  for 
Li/Si,  Jcr  would  be  higher  at  higher  Li  concentrations,  a  condition 
that  would  further  reinforce  the  conclusion  that  operation  at  higher 
concentration  is  an  effective  means  to  mitigate  failure  of  thin-film 
Li/Si  electrodes.  Such  an  analysis  with  composition-dependent  JCR, 
however,  requires  detailed  experimental  measurements  and  is  not 
attempted  in  this  paper. 

It  should  also  be  noted  that  the  analysis  here  is  based  on 
a  stress-free  initial  condition  which  does  not  consider  residue 
stresses  from  prior  lithiation  processes.  Such  a  simplification, 
however,  does  not  affect  the  primary  conclusions,  especially  the 
conclusion  that  operation  windows  at  higher  concentrations 
enhance  battery  cyclability. 

4.  Conclusions 

A  fully-coupled  finite  deformation  theory  for  analyzing  the 
coupled  mechano-diffusional  driving  forces  for  fracture  in  elec¬ 
trode  materials  is  developed.  The  formulation  of  ./-integral  treats 
the  crack  tip  as  a  separate  thermodynamic  system  and  entails 
detailed  account  of  balances  of  mass  and  energy  and  evolution  of 
entropy.  It  is  found  that  the  standard  form  of  ./-integral  for  energy 
release  rate  is  no  longer  path-independent  when  coupled 
mechano-diffusion  driving  forces  are  present.  Instead,  an  area 
integral  must  be  included  to  maintain  path-independency. 

A  numerical  scheme  for  implementing  this  path-independent 
formulation  through  mixed  finite  element  simulations  is  given. 
This  numerical  scheme  is  found  to  be  highly  robust,  even  under 
conditions  of  global  yielding. 

Under  loading,  lithium  accumulates  at  tips,  leading  to  relaxation 
of  hydrostatic  stress.  This  process  does  not  affect  in-plane  stress 
fields  or  the  tress  intensity  factor  as  along  as  the  material  is  linear 
elastic  and  sufficient  time  is  given  for  diffusion  to  occur  (i.e. 
T0S>rLi).  In  contrast,  the  deformation-diffusion  coupling  increases 
the  energy  release  rate  (driving  force  for  fracture).  This  anti¬ 
shielding  effect  on  J  can  be  understood  through  a  negative  effec¬ 
tive  Poisson’s  ratio  effect.  Specifically,  the  anti-shielding  effect 
embodies  the  fact  that  the  redistribution  of  Li  provides  another 
source  of  energy  for  the  fracture  besides  the  mechanical  fields.  Two 
complementary  mechanisms  make  this  additional  source  of  energy 
available.  These  two  mechanisms  are  controlled  by  two  partial 
atomic  volumes.  The  first  is  the  chemical-to-mechanical  flCM, 
which  controls  storage.  The  second  is  the  mechanical-to-chemical 
Q  ,  which  controls  retrieval. 

For  electrode  materials,  energy  release  rate  J  and  stress  inten¬ 
sity  factor  K  are  not  uniquely  related.  Under  conditions  of 
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significant  plasticity,  K  cannot  be  used  and  J  is  the  only  sensible 
choice  for  characterizing  fracture.  Before  global  yielding,  fracture 
is  governed  by  Irwin’s  effective  crack  length;  after  global  yielding, 
plastic  deformation  in  regions  far  from  the  crack  significantly 
reduces  J  and,  hence,  the  likelihood  of  fracture.  Global  yielding 
marks  an  inflection  point  in  the  energy  release  rate-discharge 
level  curve.  In  general,  the  driving  force  for  fracture  in  thin-film 
electrodes  can  be  lowered  by  operating  at  higher  Li  concentra¬ 
tions.  A  design  map  is  developed  to  quantify  this  effect,  accounting 
for  lithiation-induced  material  softening,  full  diffusional- 
mechanical  coupling  and  plasticity.  This  map  can  be  used  to 
guide  the  design  of  electrodes  to  improve  cyclability  and  maxi¬ 
mize  utilizable  operational  capacity. 
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Appendix  A.  Two-way  coupling  between  diffusion  and  stress 
development 


where  the  scalar  rate  factor  A  can  be  determined  from  consistency 
conditions  and  takes  the  form  of 
day 

X  =  *\D:adev  +  ^d£¥i.  (A.6) 

2<j£  2G  <Jy 

Equations  (A.5)  and  (A.6)  and  give  the  mechanical  constitutive 
behavior  when  the  yield  stress  is  reached,  while  Eq.  (A.4)  governs 
the  deformation  in  the  elastic  regime.  Together  with  the  conser¬ 
vation  of  momentum  equation 


the  above  equations  describe  the  deformation  of  the  electrode 
material. 

The  diffusive  flux  of  lithium  can  be  measured  either  in  the 
Lagrangian  frame  as  JL,or  in  the  updated  Lagrangian  frame  as  \h . 
These  measures  are  related  through  f=F-1  as 

ft  =  det(F)/^f .  (A.8) 

On  one  hand,  ]h  must  satisfy  the  requirement  of  mass  conser¬ 
vation  such  that 


To  provide  relevant  theoretical  background,  the  framework 
[15,16]  for  describing  the  constitutive  behavior  of  the  electrode 
material  is  briefly  summarized  here.  During  charge/discharge, 
lithium  atoms  diffuse  while  the  host  atoms  are  assumed  to  undergo 
only  convection  but  not  diffusion.  The  elastic  strain  associated 
with  Fe  is 


dCLi 

at 


vk[ 

dXK 


=  Rl 


(A.9) 


where  Rtf  is  the  body  source  of  Li  measured  in  the  Lagrangian 
frame.  In  this  paper,  we  assume  =  0.  On  the  other  hand,  Fick’s 
law  requires 


te— 2  ((F£’)T(Fe)-1) . 


rAii  iLi  dU  f  f  r1'8^ 

(A-1}  Jk  -  ~wMii  w 


(A.10) 


Since  s.e  is  in  general  small,  the  elastic  stretch  Ue  associated  with 
Fe  =  ReUe  is  small,  although  the  rotation  Re  may  be  finite.  The 
Cauchy  stress  <t  is  given  by 


where  D Ll  =  Du(:)  is  the  composition-dependent  diffusivity  of 
lithium,  kB  is  Boltzmann  constant;  0  is  temperature,  and  fiLi  is  the 
chemical  potential  of  lithium  which  can  be  expressed  as 


oij  =  FfK(3KEem6KL  + 206^, 


(A.2) 


where  (ft!  =  r^/3  and  rfj  =  r?.  -  are  the  isotropic  and  deviatoric 
parts  of  the  elastic  strain;  K  and  G  are  the  bulk  and  shear  modulus, 
respectively. 

Plastic  deformation  is  assumed  to  follow  the  associated  flow 
rule  with  a  composition-dependent  yield  surface  ay  =  oY(CLi)  in 
the  form  of 

0.  =  Q(o)  =  l<7dev  ■■  odev  -  \(oY)2  =  0,  (A.3) 

where  adev  as  -  tr(a) /3I  is  the  deviatoric  part  of  the  Cauchy  stress. 
Below  the  yielding  threshold,  D  [Eq.  (2)]  is  simply  given  by  the  rate 
form  of  Eqs.  (3)  and  (A.2),  i.e. 


"“-ntl  ca.1i) 

Here,  am=akk/3  is  the  hydrostatic  stress,  QLl(SF^  =  dJSF /dCu  is  the 
stress-free  partial  atomic  volume  of  lithium,  and  0  =  0(Fe,Cil,0)  is 
the  Lagrangian  density  of  the  Helmholtz  free  energy.  In  this  paper, 
we  consider  the  free  energy  to  be  of  the  form 


=  <t>SF(cLi„d) 

w^K(  4)2+G, 


(A.  12) 


In  the  above  relations,  ij>SF  and  w  are,  respectively,  contributions 
to  0  due  to  stress-free  chemical  interactions  and  mechanical 
stresses.  Under  the  assumption  that  elastic  strain  re  is  small,  Eq. 
(A.11 )  can  be  further  simplified  into 


D  = 


tr(“")1 


=  n^.  -  det(Fe)f/1(SF,<7m 

(A.4)  Here, 


(A.13) 


Here,  W-cr+tT-W  is  the  objective  Jaumann  rate  of  <t, 

with  W  =  (L-Lr) /2  being  the  spin  tensor.  Upon  yielding,  the  rate  of 
deformation  including  plasticity  is 


P-sf  =  WSF/dCu  =  /4(CU)  (A.14) 

is  the  isothermal  chemical  potential  of  lithium  at  zero  stress. 
Generally  speaking,  n^F(CLl)  can  be  determined  by  fitting  the 
activity  coefficient  to  experimentally  measured  open  circuit 
potential  (OCP)  data  [2],  To  simplify  the  problem,  however,  we 


m  /  Journal  of  Power  Sources  230  (2013)  176-193 


adopt  the  ideal-solution  form  of  fi{fF  originally  proposed  by  Larche 
and  Cahn  [60]  and  subsequently  used  by  Purkayastha  and 
McMeeking  [61]  and  Yang  et  al.  [62]  in  the  form  of 


/4f  =  F*  +  Min 


?max  -  V 


(A.15) 


where  ?max  is  the  maximum  charge  limit  of  the  electrode 
material  and  constant  $  is  a  reference  chemical  potential.  The 
composition-dependent  diffusivity  DLi,  on  the  other  hand,  is  taken 
to  be  [61] 


DLi  =  Dq(1  -  2). 


Therefore,  invoking  the  divergence  theorem  ( 


Q 

The  Lee  decomposition  [Eq.  (1)]  indicates  that 


y  = 


dFf! 


(B.4) 


(A.16) 


ainVSF 

9X|  “ax/  'if  1  ‘ir  9X|  ‘ '*»  0X1  ’ 

For  free  energy  of  the  form  0  =  0(F e,Ci!,0),  the  1st 
Piola-Kirchoff  stress  and  Li  chemical  potential  are,  respectively, 


f  W  -  $\  and 

I  0tjk \c“,0 

I  aLi  =  _  ^t\  Fe9lnySF| 

[ M  IPV,Fj9  dFl\cufi  *  sc«  U  ' 


-ill  f-0«|81nVgf. 

'0C«|Pe  IC«  ’ 


Here,  S=?/fmax  =  CLi/C^3X  is  the  state-of-charge  (SOC),  which 
represents  the  amount  of  Li  relative  to  the  fully-charged  state. 
Combining  Eqs.  (A.10— A.16)  leads  to  the  following  form  of  the 
lithium  flux, 


wm  +  ttmech ,  with 

ar Li 

rtjfKiffiM]  and 


M  M]t  ax/ 


In  the  above  relations,  J^'chem  and  JLpmech  are  contributions  due 
to  chemical  interaction  and  mechanical  driving  force,  respectively. 
According  to  Purkayastha  and  McMeeking  [61  ],  the  main  advantage 
of  Eqs.  (A.15)  and  (A.16)  is  that  it  correctly  captures  the  fact  that  the 
stress-gradient-driven  diffusion  flux  J^'mech  [cf.  Eq.  (A.17)]  vanishes 
at  the  charging  limits  of  ?  =  0  and  c  =  ?max. 

A  mixed  finite  element  framework  [29,30]  is  used  to 
analyze  the  coupled  diffusion/deformation  processes  delineated 
by  Eqs.  (A.3)— (A.17).  Details  of  this  numerical  framework  are  given 
in  [29], 

Appendix  B.  Proof  of  equation  (24) 

J*(T)  as  defined  in  Eq.  (23)  is  path-dependent  and  equals  to 

J*m  =  ~  J wpSyNjdr0  +  J ^'dAo  (B.l ) 


for  any  continuum  domain  Q  bound  by  boundary  T  =  a Q.  To 
demonstrate  this,  we  first  note  that  balance  of  momentum  [Eq. 
(A.7)]  requires  that 


where  fp  =  (Fp)  1 .  Therefore,  under  isothermal  conditions, 

00  00  I  0Fy  00  I  0CLi 

0X,  =  0Fg|cu0Xi+0C«|raXi 

L  vsfFp  <tpk1  dF*j  +  (uLi  +  F  (TPK1 9ln'/SFA  dCb  (B.6) 

-  r  0Xt  +  x  +FjK<Tk)  0C«  )  0Xj 

=  vSFFPaPKidFV  |  uu»Cu  p  pK]0lnl/Sf 

v  axt  +  M  ax,  +  F]K  kj  ax. 

In  the  derivation  of  Eq.  (B.6),  it  is  assumed  that  VSF  is  uniquely 
determined  by  current  Li  concentration  [Eq.  (3)]. 

Combining  Eqs.  (B.4)  and  (B.6)  leads  to 

(BJ) 

With  the  assumption  that  elastic  strains  are  small  (i.e.,  FjjF! j  =  <5J(), 
the  second  term  on  the  right  hand  side  of  Eq.  (B.7)  can  be  easily 
transformed  into 

(B-81 

where  the  plastic  strain  is  defined  as  As-^FJ^  -  8jK).  By  Eqs. 
(B.3),  (B.7)  and  (B.8),  one  obtains  2 

J*m  =  j /'^dAo  -  J  (vSF)2app^dAQ.  (B.9) 

Q  Q 

To  proceed  further,  we  define  the  plastic  work  as 

wp  =  J  det(F)a  :  I fdt~  J  (vSF)2 <r$2EjKdt 


(B.10) 
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and  adopt  the  standard  approximation  such  that 


Substitution  of  Eq.  (B.ll)  into  Eq.  (B.9)  and  the  divergence 
theorem  lead  to 


/( 90)  =  J ~  J wpdyNjdr0.  (B.12) 

a  a  Q 

The  validity  of  Eq.  (B.12)  mostly  depends  on  how  significant 
the  error  associated  with  approximation  (B.ll)  is.  The  path 
independency  of  Eq.  (30)  is  therefore  less  reliable  when  approxi¬ 
mation  (B.ll)  breaks  down.  Indeed,  for  cracks  in  elasto-plastic 
materials,  contours  that  are  very  close  to  the  crack  tip  should  not 
be  used. 


Appendix  C.  Implementation  of  J  using  the  energy  domain 
integral  method 

Equation  (30)  is  the  path-independent  form  of  the  energy 
release  rate  used  in  this  paper  for  the  numerical  evaluation  of  G.  To 
facilitate  numerical  calculation,  note  that  Eq.  (A.ll)— (A.13)  allow 
Eq.  (26)  to  be  further  simplified  into 

J  =  m  =  j[(w  +  wp)Sv-  afP  (dUi/dX, )]  Njdr0 

“/( ~  det(r)QL,(SF>am)  fx,  dA°- 

which  consists  of  a  contour  integral  over  a  finite  path  and  an  area 
integral  over  a  finite  domain. 

Based  on  Eq.  (C.l ),  we  adopt  the  method  originally  proposed  by 
Shih  et  al.  [40]  in  order  to  facilitate  the  computation  of  J(T).  For 
arbitrary  contours  r  and  C0  (cf.  Fig.  1),  consider  a  sufficiently 
smooth  field  of  q*  =  q*(X)  such  that  q  =  1  on  T  and  q*  —  0  on  Co- 
Here,  f  is  not  necessarily  infinitesimal.  For  traction-free  crack 
surfaces  C_  and  C+,  the  path-independent  expression  of  energy 
release  rate  [Eq.  (C.l )]  can  be  rewritten  as 

J  =  Kr )  =  J  [~  (w  +  Wp) 5y  +  oj^1  (Otii/aXi)]  q*Njdr o 

(C.2) 

Here,  —T  stands  for  contour  T  with  the  opposite  loop  direction. 
Application  of  the  divergence  theorem  leads  to 

1-m  -  fr[-a(wdx,w')+»r^]‘»° 

+j%[-('“+wry'i+atM-]dA°  <C3) 


With  Eqs.  (A.ll)— (A.13)  and  (C.3), 

,  +| -  *1 J  (C.4) 

+/,-[detd»)C“lSFV„-g|F  J  gtv 

Here,  T  and  C0  are  arbitrary  (cf.  Fig.  la)  and  F  is  not  necessarily 
infinitesimal. 


Appendix  D.  The  path  independency  of  J  and  associated 
numerical  error 


The  accuracy  of  the  path-independent  expressions  of  J  [Eqs.  (30) 
and  (31 )]  is  contingent  upon  the  accuracy  of  the  approximation  in 
Eq.  (B.ll).  To  validate  this  approximation,  we  consider  two  inte¬ 
gration  paths  as  shown  in  Fig.  D.l(b).  The  material  parameters  and 
system  dimensions  are  the  same  as  those  used  for  Fig.  9,  with 
oY  =  0.01. 


Up  to  the  point  of  ASy,  fp(Ffar)  and  jcp(rnear )  [and  also 
jNCP(  ffar)  ancj  jNCP^near-jj  caicuiated  for  the  two  paths  shown  in 
Fig.  D.l(b)  differ  by  less  than  0.5%.  However,  this  error  begins  to 
increase  once  the  film  yields  globally.  This  divergence  is  intrinsic  to 
all  calculations  of  ./-integral  and  is  not  specific  to  the  diffusion- 
deformation  two-way  coupled  problems  considered  here.  Even 
for  purely  mechanical  cases,  J  calculated  for  different  paths  do  not 
show  good  asymptotic  behavior  when  r~*  0.  Instead,  /  converges  to 
the  same  value  of  J“  when  r->°°[63].  Therefore,  although  by 
definition  per  se  the  energy  release  rate  is  given  with  T->0 
[Eq.  (20)],  computations  should  be  performed  on  T  that  is  as  far  as 
possible  from  the  crack  tip.  We  therefore  perform  all  our  calcula¬ 
tions  for  /  using  the  largest  possible  contour.  For  all  calculations 
involved  in  this  paper, JCP(rfar)  and  fP(fnear)  never  differ  by  more 
than  10%. 
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Nomenclature 


a:  Crack  length  [m] 

cH,  cLi:  Eulerian  concentration  of  host  (H)  and  lithium  (Li)  [atoms  m-3] 

CH,  C,J :  Lagrangian  concentration  of  host  (H)  and  lithium  (Li)  [atoms  m-3] 
rfcmck.  Nominal  out-of-plane  thickness  of  the  2-D  cracked  system  [m] 

D,  De,  DSF,  Dp:  Total,  elastic,  stress-free  volumetric  and  plastic  rates  of  deformation 

.  I5-’] 

DLl:  Concentration-dependent  diffusivity  of  Li  [m2  s  1  ] 

D'q1:  Characteristic  concentration-independent  diffusivity  of  Li  [m2  s-1] 
e:  Internal  energy  per  unit  reference  volume  [J  m-3] 

E:  Young's  modulus  [Pa] 

F,  F,  Fsf,  Fp:  Total,  elastic,  stress-free  volumetric  and  plastic  deformation 
gradients 

G:  Shear  modulus  [Pa] 

J:  Energy  release  rate  for  fracture  [J  m-2] 

Jcr:  Critical  energy  release  rate  for  fracture  [J  m-2] 

JSF:  Volume  change  due  to  changes  in  Li  concentration  [/SF  =  det(FSF)[ 

J11:  Lagrangian  flux  of  lithium  [atoms  m-2  s_1] 

Jq:  Lagrangian  flux  of  heat  [J  m-2  s-1] 

J’:  Lagrangian  flux  of  entropy  [J  K-1  m-2  s-1] 

K:  Bulk  modulus  [Pa] 


K,  :  N 


jr [Pa  m 


m'/2] 


re  toughness  [Pa  m 

l:  Crack  length  measured  in  Lagrangian  configuration  [m] 
nu:  Amount  of  Li  in  the  crack  system  per  unit  d€rack  [atoms  m-1] 

Scrack:  Entropy  of  crack  subsystem  per  unit  dFack  [J  K_1  m-1] 

To:  Characteristic  loading  time  [s] 

Ucrack:  Internal  energy  of  crack  subsystem  per  unit  dcrack  [J  m”1] 
uu:  Partial  internal  energy  of  Li  [J  per  atom] 
v:  Velocity  of  material  point  [ms-1] 

VSF:  Isotropic  stretch  due  to  changes  in  Li  concentration  [VSF  =  (/SF)1/3[ 
w:  Elastic  energy  per  unit  Lagrangian  volume  [J  rrr3] 
wp:  Plastic  potential  [J  m-3] 

x;  Position  of  material  point  in  current  configuration  [m] 

X:  Position  of  material  point  in  reference  configuration  [m] 


Surface  energy  per  unit  crack  area  [J  m-2] 

Ee:  Elastic  strain 

feg:  Effective  plastic  strain 

n:  Entropy  per  unit  reference  volume  [J  K-1  m-3] 

8:  Temperature  in  the  continuum  [K] 
gcrack.  xemperature  in  the  crack  subsystem  [K] 

HLi:  Chemical  potential  of  Li  in  the  continuum  [J  per  atom] 
jjLi,crack.  chemical  potential  of  Li  in  the  crack  subsystem  [J  per  atom] 
n'j:  Reference  chemical  potential  of  Li  [J  per  atom] 
v :  Poisson's  ratio 

l:  Lithium  number  per  host  in  chemical  formula  LijH 
fmax-'  Maximum  5  at  the  fully  charged  state 
2:  State  of  charge  (SOC)  [2=f/fmax] 

20:  Initial  SOC 
2:J3omain-average  SOC 

ASy:  Discharge  level  at  which  thin  film  begins  to  yield  globally 
AScr:  Critical  discharge  level  at  which  j  =  Jcr 
a:  Cauchy  stress  [Pa] 

<rm,  <rmises:  Hydrostatic  and  von  Mises  invariants  of  Cauchy  stress  [Pa] 
<Tpin,  aPK2:  First  and  second  Piola-Kirchhoff  stresses  [Pa] 


